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1  Introduction 


One  of  the  fundamental  and  powerful  ideas  of  signal  processing  is  that  of  a  system  function  and 
input-output  relations.  This  is  commonly  symbolized  in  the  following  way.  If  we  have  an  input  time 
function,  f(t)  which  passes  through  a  system  characterized  by  a  system  function  h(t),  then  the  output 
is  given  by  x(t).  This  is  symbolized  in  Fig.  1. 


At) 


h(t) 


x(t) 


Figure  1:  Input-output  representation  of  a  system  in  the  time  domain. 


Equivalently  in  the  Fourier  domain,  if  the  input,  system  and  output  transforms  of  the  time  functions 
are  given  by  F(uj),  H(u)  and  X(oj)  respectively,  then  the  input-output  relations  are  symbolized  as  in 
Fig.  2. 


F(co) 


H(co) 


X(co) 


Figure  2:  Input-output  representation  of  a  system  in  the  frequency  domain. 


However,  over  the  last  fifty  years  it  has  been  found  that  many  natural  and  man  made  signals  are 
nonstationary,  and  the  above  formulation  does  not  fully  describe  what  is  happening.  Our  fundamental 
point  to  be  developed  in  this  report  is  that  if  the  signals  are  nonstationary  an  immense  simplification 
and  advantage  occurs  if  we  have  input-output  relations  in  the  time-frequency  plane.  We  symbolize 
this  in  Fig.  3.  where  Cf(t,uj)  and  Cx(t,  ui)  (e.g.  spectrogram,  Wigner,  etc.  distribution)  are  the  input 
and  output  time-frequency  distributions  and  the  box  with  the  question  mark  is  meant  to  symbolize  the 
time- frequency  system  function.  It  is  one  of  the  aims  of  this  report  to  explain  how  the  time-frequency 
system  function  can  be  obtained  and  to  show  the  advantages  of  such  a  formulation. 

To  illustrate  and  motivate  our  method  we  start  with  a  simple  example.  Consider  the  differential 
equation 

+  =  (1) 
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Cf{t,co) 


Cx  ( Uo ) 


Figure  3:  Input-output  representation  of  a  system  in  the  time-frequency  domain. 


where  f(t)  is  a  given  driving  force.  Perhaps  there  is  no  more  studied  equation  than  this  one.  In 
principle  this  equation  can  be  solved  “exactly”  by  many  methods.  For  both  practical  reasons  and  to 
gain  insight,  one  often  transforms  this  equation  into  the  Fourier  domain.  Defining 


X(u)  = 

F(v)  = 


1 


\/2vr 


Sx(t)^dt 


f(t )  e~ltuJ  dt 


the  differential  equation  transforms  into 

[—uj2  +  2  inu  +  ujq\X(u>)  =  F(u) 

whose  exact  solution  is 


X{u)  = 


F(w) 


(2) 

(3) 

(4) 

(5) 


[— uj2  +  2  i/Mjj  +  0Jq] 

The  reasons  for  going  into  the  Fourier  domain  are  many.  First,  it  offers  a  practical  way  of  solution 
since  now  one  can  find  the  time  solution  by  way  of 


x(t)  = 


1 


n») 


situ  dt 


(6) 


a/27 t  J  [—(X2  +  2 ifj,LU  +  Wq] 

Perhaps  more  importantly  the  reason  for  going  into  the  Fourier  domain  is  that  one  can  gain  insight  into 
the  nature  of  the  solution  and  both  reasons  have  become  part  of  standard  analysis  in  both  engineering 
and  physics,  as  exemplified  by  input-output  relations. 

Now  consider  a  specific  example.  Take 


d2x(t)  n  dx(t)  2  -at2 /2+i0t2 /2+iuiit 

dt2  +  ^  dt  +UJ°X~e 


(7) 


The  solution  obtained  numerically  is  shown  in  Fig.  4.  The  Fourier  transform  of  the  driving  force  is 


F  M  = 


1 


\Ja  —  i/3 


exp 


a(uj  —  uj\)2  ,P(uj  —  LOi)- 


2(a2  +  (32)  2 (a2  +  (32)  _ 


which  gives 


X{u)  = 


exp 


i)2  ■ uji)2 

'  2(o2+/32)  Z2 (a2+/32) 


y/a  —  i/3  [ —uj 2  +  +  Wq] 


(8) 


(9) 
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x  10" 


t(s) 


Figure  4:  Solution  of  Eq.  (351),  x(t).  The  parameters  are  /i  =  1,  w0  =  67t  rad/as,  a  =  .001,  (3  =  6/57T, 
u>\  =  —  87t  rad/as. 


Figure  5:  Frequency  spectrum  of  x(t)  shown  in  Fig.  4.  The  two  peaks  are  due  to  the  resonances  of  the  oscillator 
located  at  /  =  ±3  Hz. 

and  the  squared  magnitude  of  this  is  shown  in  Fig.  5. 

In  Figs.  4-5  we  plot  the  signal  and  spectrum  for  the  values  indicated  in  the  caption.  Much  can 
be  learned  from  a  study  of  x{t )  and  X(u>).  However,  even  more  can  be  learned  than  is  commonly 
discussed  in  textbooks  as  we  now  show  if  we  plot  the  time-frequency  distribution.  In  Fig.  6  we  plot 
a  possible  C(t,u)  for  the  signal  x(t).  We  see  that  something  remarkable  happens:  one  gets  a  simple, 
clear  picture  of  what  is  going  on  and  of  the  regions  which  are  important  [?].  Such  distributions  have 
been  studied  for  over  seventy  years  in  the  field  of  time- frequency  analysis  in  engineering  [?,  ?,  ?],  but 
the  system  function  approach  has  not  been  developed 

In  this  report  we  consider  systems  that  are  described  by  differential  equations.  For  the  sake  of 
clarity  we  will  developed  the  ideas  for  ordinary  differential  equations  with  constant  coefficients  and 
then  give  the  other  cases  in  the  Appendix.  However  in  the  next  Section  we  first  give  a  brief  review  of 
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123456789  10 

t(s) 


Figure  6:  Time-frequency  distribution  of  x(t)  represented  in  Fig.  4.  The  main  energy  response  occurs  when 
the  forcing  function  hits  the  resonant  frequency  of  the  oscillator,  which  is  located  at  /  =  3  Hz.  Note  that  the 
picture  shows  only  positive  time  and  frequencies. 

t ime- frequency  dist ribut ions . 

2  Brief  Review  of  Time-frequency  distributions 

Starting  in  the  early  1940s  it  was  realized  that  for  many  natural  and  man  made  signals  their  spectra 
change  over  time.  The  development  of  the  physical  and  mathematical  ideas  needed  to  explain  and 
understand  time-varying  spectra  has  evolved  into  the  field  now  called  “time-frequency  analysis”  or 
time-varying  spectral  analysis.  The  purpose  of  time-frequency  analysis  is  to  understand  the  nature  of 
signals  so  that  we  may  understand  the  physical  phenomena  generating  them,  the  medium  of  propaga¬ 
tion,  their  structure,  classification,  detection,  etc.  The  basic  idea  is  to  find  a  joint  density  of  time  and 
frequency  that  indicates  what  frequencies  are  present  in  the  signal  and  how  they  are  changing  in  time. 
The  main  initial  impetus  occurred  at  Bell  Laboratories,  with  the  development  of  the  “sound  spectro¬ 
graph”  for  studying  speech.  Sometime  later  it  was  realized  that  the  principle  methods  used  in  engi¬ 
neering  for  time-frequency  analysis  are  mathematically  analogous  to  the  development  that  occurred 
in  quantum  mechanics  starting  with  Wigner,  and  now  called  the  phase  space  formulation  of  quantum 
mechanics  independently  of  their  use  in  quantum  mechanics,  these  distributions  have  been  developed 
in  signal  analysis  as  a  means  of  understanding  how  the  spectral  content  of  a  signal  changes  in  time  for 
classical  variables.  This  development  has  a  long  history  and  originated  with  the  work  of  Koenig,  Dunn, 
and  L.  Lacy  XciteKoenig,  Gabor  XciteGabor,  and  Ville  XciteVille,  and  since  that  time  there  have 
been  many  quasi-distributions  or  representations  that  have  been  used  and  developed.  A  general  for¬ 
mulation  of  quasi-distributions  was  given  by  Cohen  Xciteleon3,  and  many  methods  have  been  devised 
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for  obtaining  distributions  with  desirable  properties  Xcitechoi,Zhao,Jeong,Salal,Sala2,Amin2,pat3. 

Of  particular  relevance  to  our  considerations  is  the  wide  range  of  applications  of  these  distribu¬ 
tions  to  classical  systems,  that  includes  acoustics  XciteGaun2,pat4,  speech  processing  XcitePitton, 
musical  instruments  XcitePielemeier,  machine  monitoring  XciteRizzoni, Atlas,  stochastic  processes 
XciteAmin3,pitton2,  biomedical  signals  XciteWilliams,Wood,  sonar  and  radar  XciteGaunl,  nonlin¬ 
ear  dynamic  al  systems  XciteLorenzo,  among  many  others  Xcitepat.  The  fundamental  idea  of  this 
approach  is  to  study  the  time-frequency  properties  of  a  classical  variable  such  as  a  pressure  wave, 
current,  etc.  The  distributions  are  typically  calculated  from  experimental  data  or  by  first  generating 
the  numerical  solution  from  the  governing  differential  equation. 

2.1  Mathematical  Development 

The  development  of  quasi-distributions  or  time-frequency  representations  occurred  more  or  less  simul¬ 
taneously  in  both  quantum  mechanics  and  signal  analysis  although  from  very  different  perspectives 
[?,  ?].  The  basic  objective  is  to  devise  a  joint  density  in  time  and  frequency.  One  can  set  up  the  issue 
as  follows.  If  we  have  a  signal  s(t)  and  its  Fourier  transform  S(uj),  then  the  instantaneous  power  is 

|  s(t)  | 2  =  intensity  per  unit  time  at  t 
and  the  density  in  frequency,  the  energy  density  spectrum,  is 

|  S(uj)  | 2  =  intensity  per  unit  frequency  at  uj 
What  one  seeks  is  a  joint  density,  P(t,uj),  so  that 

P(t,  uj)  =  the  density  (or  intensity)  at  time  t  and  frequency  uj. 


Ideally  the  joint  density  should  satisfy 


J  P(t,  uj)  du> 


j  P(t,  uj)  dt 


s(t)  I2 

(10) 

S(uj)  |2 

(11) 

which  are  called  the  time  and  frequency  marginal  conditions.  Wigner,  and  later  Ville,  gave  such  a 
function,  which  is  now  called  the  Wigner-Ville  distribution  [?,  ?], 

W(t,uj)  =  —  f  s*(t  —  t/2)  s(t  +  r/2)  e_?'™  dr  (12) 

27t  J 

It  satisfies  the  time-frequency  marginal,  Eqs.  (10)  and  (11),  but  in  addition  it  has  the  property  that 
the  first  conditional  moment  of  time  at  a  given  frequency  is  given  by  the  instantaneous  frequency  of 
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the  signal.  An  important  associated  concept  is  the  instantaneous  frequency  requirement.  If  we  write 
a  signal  in  terms  of  amplitude  and  phase  as 


s(t)  =  A(t) 


(13) 


one  takes  the  instantaneous  frequency  to  be 


UJi{t)  = 


dip(t) 

dt 


Now,  the  conditional  frequency  for  a  given  time  is  given  by 

1 


= 


s(t) 


LO  P(t,  uj)  du ) 


For  many  distributions,  such  as  the  Wigner  distribution 


{i0)t=  UJi(t)  = 


dcp(t) 

dt 


(14) 


(15) 


(16) 


From  this  point  of  view  the  instantaneous  frequency  can  be  thought  of  as  the  average  frequency  at 
a  particular  time.  We  point  out  that  for  certain  distributions,  Eq.  (16)  is  exactly  satisfied  and  also, 
often,  the  peak  of  the  distribution  is  approximately  given  by  .  It  is  important  to  mention  though 
that  writing  a  real  signal  in  a  complex  form  as  given  by  Eq.  (13)  has  a  long  history.  The  basic  issue  is 
that  there  are  an  infinite  number  of  ways  to  generate  a  complex  signal  from  a  real  signal.  Currently 
the  most  widely  accepted  method  is  the  one  devised  by  Gabor  and  called  the  analytic  signal.  We  do 
not  address  these  issues  here.  This  is  not  true  for  all  distributions,  but  is  so  for  the  Wigner  and  others. 


Among  other  quasi-distributions  subsequently  proposed  in  signal  analysis  and  quantum  mechanics 
were  the  Rihaczek,  Page,  and  Margenau-Hill,  among  others.  In  1966  a  method  was  devised  that  could 
generate  in  a  simple  manner  an  infinite  number  of  new  ones  [?].  This  general  class  is  given  by 

C(t,u)  =  -^2  JJjs*{u-T/2)s(u  +  T/2)(j>(0,T)e-m-iT“+ieududTde  (17) 

where  </>($,  r)  is  a  two  dimensional  function  called  the  kernel.  The  properties  of  a  distribution  are 
reflected  as  simple  constraints  on  the  kernel,  and  by  examining  the  kernel  one  readily  can  ascertain 
the  properties  of  the  distribution.  This  allows  one  to  pick  and  choose  those  kernels  that  produce 
distributions  with  prescribed,  desirable  properties.  Williams  and  co-workers  devised  and  crystallized 
the  idea  of  kernel  design.  They  developed  a  methodology  for  the  construction  of  densities  with 
desirable  properties.  At  about  the  same  time,  Zhao,  Atlas  and  Marks  similarly  produced  a  density 
that  resolved  many  of  the  difficulties  with  the  Wigner  distribution.  These  works  and  others  led  to 
major  developments  in  trying  to  understand  the  nature  of  these  time-frequency  densities  and  also  to 
practical  applications  to  the  fields  mentioned  above  [?,  ?,  ?,  ?]. 
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Figure  7:  The  Wigner  distribution  of  s(t)  =  e  at2 /2+lP^ /'2+lu}ot  shown  in  the  central  part.  The  top 

figure  is  the  real  part  of  the  signal  and  the  left  figure  is  the  absolute  of  the  spectrum. 

One  can  think  of  a  time-frequency  distribution  as  a  two  dimensional  transform  of  a  one  dimensional 
function.  In  the  mathematical  sense  the  distribution  contains  the  same  information  as  the  signal  since 
it  is  constructed  from  it  and  for  many  distributions  the  signal  can  be  obtained  from  it  uniquely. 
However  a  dramatic  thing  happens  when  one  plots  or  studies  the  time-frequency  distribution  instead 
of  the  signal  or  spectrum,  namely  the  physical  nature  of  the  signal  becomes  much  clearer.  The  best 
way  to  understand  the  concept  of  a  time- frequency  representation  is  to  consider  a  number  of  examples. 

2.2  Example  1 

Consider  the  signal, 

s(t)  =  e~at2 /z+ifit2 /2+iuJo  t  (18) 

In  Fig.  7  we  plot  the  distribution.  In  the  top  panel  is  the  real  part  of  the  signal,  the  left  panel  is  the 
absolute  value  of  the  spectrum  and  in  the  main  figure  we  have  the  time-frequency  distribution,  and 
in  this  we  show  the  Wigner  distribution.  Notice  that  it  is  totally  concentrated  along  the  curve 

oj  =  u>o  +  (3t  (19) 

which  is  exactly  the  instantaneous  frequency. 


2.3  Example:  Multipart  signals 

One  of  the  interesting  aspects  of  time-frequency  analysis  is  that  it  reveals  when  a  signal  consists  of 
parts.  We  use  the  word  parts  instead  of  components  since  sometime  “components”  is  often  associated 
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with  orthogonal  components  which  is  not  necessarily  the  case  here.  Consider  the  signal 
=  g-«l*2/2+^2 1  _j_  e— a\t2 /2+if)\t2 /2+iu\t.  _|_  at2 /2+i^t3 /3+if3t2 /2+iuj0t 


(20) 


In  Fig.  8  we  show  the  Wigner  distribution  of  s(t).  Note  how  clearly  the  time-frequency  distribution 
reveals  that  it  consists  of  three  parts  and  moreover  note  that  each  part  is  in  some  sense  approximately 
concentrated  along  its  instantaneous  frequency,  respectively  uq(f),  uj2(t)  and  0J3 (t) 


U>i(t)  =  i02 

(21) 

W2(t)  =  Ul  +  fat 

(22) 

(t)  =  'yt2  +  f3t  +  uj0 

(23) 

Example:  whale  sound 

For  a  first  example  we  take  a  whale  sound.  In  Fig.  9,  running  across  the  page,  above  the  main 
figure  is  the  sound,  that  is  the  air  pressure  as  a  function  time.  To  the  left  of  the  main  figure  is  the 
energy  density  spectrum,  that  is,  the  absolute  square  of  the  Fourier  transform  of  the  signal.  The 
energy  density  spectrum  indicates  what  frequencies  existed  and  what  their  relative  strengths  were  for 
the  duration  of  the  signal  but  it  does  not  indicate  when  these  frequencies  occurred.  For  this  sound 
it  tells  us  that  the  frequencies  ranged  from  about  175  to  about  375  cycles  per  second.  The  main 
figure  is  a  time-frequency  plot  and  now  we  can  determine  not  only  what  frequencies  occurred  by  when 
they  occurred  and  what  were  their  relative  intensities  were  as  time  evolves.  At  the  start  the  starting 
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Figure  9:  A  whale  sound. 


frequency  was  about  175  Hz  and  increased  more  or  less  linearly  to  about  375  Hz  in  about  half  a  second, 
stayed  there  for  about  a  tenth  of  a  second,  and  then  decreased  approximately  linearly  to  about  200 
Hz. 

Now  consider  Fig.  10  where  the  signal  is  made  up  of  four  sine  waves  with  varying  duration.  If  one 
were  to  look  at  the  spectrum  the  only  conclusions  one  could  arrive  at  is  that  four  constant  frequency 
signals  existed,  but  with  no  information  about  when  they  existed.  The  time-frequency  plot  does  that 
so  clearly  and  also  shows  the  amplitude  variations. 

Now  consider  the  signal 

s(t)  =  e-Qii2/2-M wit  _|_  e~a2t2 /2+i(3t2 /2+iuj2t  (24) 


The  distribution  is  plotted  in  Fig.  11.  It  is  concentrated  along  the  instantaneous  frequency  of  each 
part 

u  =  ujq  +  /3t  ;  uj  =  loq  (25) 

We  call  such  signals  Multipart  because  they  consists  of  parts,  in  this  case  two.  One  of  the  advantages 
of  time-frequency  analysis  is  that  it  effectively  shows  that  a  signal  consists  of  parts,  something  that 
can  not  be  seen  directly  in  the  signal  or  spectrum. 

As  a  last  example  consider 

s(t)  =  (ck/tt)1/4  e-at2/2+^4/4+nt3/3+^t2/2+ia;oi  (26) 
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Figure  10:  Four  constant  frequency  sine  waves  of  different  duration. 
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Figure  12:  The  Wigner  distribution  of  Eq.  (26). 

whose  distribution  is  illustrated  in  Fig.  12.  Again  the  concentration  is  along  the  instantaneous 
frequency 

ui  =  r]t3  +  7 f2  +  (3t  +  u>  o  (27) 


3  Outline  of  the  Approach 

Our  aim  is  to  develop  input-output  relations  for  nonstationary  situations  and  when  the  governing 
equations  relating  input  and  output  is  a  differential  equation.  Although  the  details  will  change  when 
dealing  with  ordinary  and  partial  differential  equations,  the  general  approach  is  conceptually  the 
same,  and  hence  we  outline  the  method  in  broad  terms.  Details  can  be  found  in  the  published  papers. 
Suppose  the  governing  equation  is 

L[u\  =  f  (28) 

where  L  is  a  linear  operator  and  /  is  a  known  driving  force.  This  situation  covers  both  the  ordinary 
and  partial  differential  equation  possibility.  In  the  former  case  u  is  x(t),  and  in  the  case  of  a  field  u 
would  be  a  function  of  x  and  t.  that  is  u(x,t). 

We  will  now  show  the  approach  in  the  case  of  ordinary  differential  equations.  The  aim  is  to  obtain 
a  differential  equation  for  the  Wigner  distribution  of  x[t)  which  will  involve  the  Wigner  distribution 
of  /(f),  and  therefore  the  starting  equation  is 

L\x(t)\  =  /(f)  (29) 
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Step  1:  Take  the  Wigner  distribution  of  both  sides  to  obtain 

WL[x\,L[x](t,u)  =  Wfj(t,u)  (30) 

We  now  have  to  express  W£[xi  Mxi(t,  lo)  in  terms  of  Wx>x(t,  u). 

Step  2:  As  is  often  the  case  L  itself  is  the  sum  of  operators, 


II 

3 

(31) 

giving 

^2WLn[X],Lm[x}(t,^)  =  Wfj(t,Uj) 

n,m 

(32) 

Step  3:  Express 

WLn[xiLrn[x]{t,uj)  in  terms  of  Wx>x{t,u) 

(33) 

This  will  be  discussed  in  detail  for  a  variety  of  operators. 

Other  Equations.  It  is  also  of  some  interest  to  consider  equations  of  motion  for  the  cross  Wigner 
distributions.  Starting  with  Eq.  29  write 

L[x{t+  \t)\  =  f(t+  \t)  (34) 

and  multiply  both  sides  by  x*(t  —  ^r)e_*TaJ  and  integrate  with  respect  to  r.  We  have  therefore 

J  x*  (t  —  \t)  L[x(t  +  \t)\  e~tTUJ dr  =  j  x*  (t  —  \r)  f  (t  +  ^r)  e~lTU dr  (35) 

which  gives 

WXtL[x](t,u)  =  Wxj(t,u)  (36) 

Similarly, 

WL[  x],x(t,u)  =  Wf}X(t,iv)  (37) 

These  equations  can  be  used  to  manipulate  intermediate  results  in  the  derivation  of  the  equation  for 

Wx,x. 

4  Transformation  of  Ordinary  Differential  Equations  with  constant 
coefficients  into  phase  space 

Many  input-output  linear  systems  are  characterized  by  ordinary  differential  equations  and  the  study  of 
such  systems  is  fundamental  in  a  number  of  branches  in  signal  processing  [?] .  The  differential  equation 
may  or  may  not  have  time-dependent  coefficients.  It  is  often  the  case  that  these  systems  exhibit  time- 
varying  frequencies,  but  generally  speaking  time- frequency  methods  have  not  been  applied  directly  to 
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these  types  of  systems.  It  has  been  done  only  in  a  circuitous  way,  namely  by  first  numerically  solving 
the  differential  equation  and  then  substituting  the  solution  into  a  time-frequency  distribution  [?].  For 
the  input  signal  we  use  f(t)  and  for  the  output  x(t),  while  for  the  governing  equation  we  take 


dnx  dn  lx  dx 

■■■  +  a1—  +  a0x  =  f(t) 


or,  in  polynomial  notation 


p(D)X(t)  =  m 


where 

P(D)  =  anDn  +  an_iDn_1...  +  aiD  +  a0 


The  Wigner  distribution  is  defined  by  XciteWigner 


(38) 


(39) 

(40) 


WXjX{t,u)  =  -^  J  x* (t  -  \t)  x(t  +  \t)  e  lTU}  dr  (41) 

Now,  take  the  Wigner  distribution  of  both  sides  of  Eq.  (39)  to  obtain 

W P(D)x,P(D)x{t->UJ )  =  (42) 

We  now  state  our  main  result:  the  input-output  governing  equation  for  the  Wigner  distribution 
associated  with  Eq.  (38)  is  given  by 

P*(A)P(B)WXjX(t,u)  =  Wfj(t,Lo)  (43) 

where 

Id  Id 

A=~~ju  B  =  -i-+ju  (44) 

2  dt  J  2  dt  J  v  ' 

and  the  star  sign  stands  for  complex  conjugation  of  the  constants  ao,--.,an.  The  distribution 

Wfj(t,uj)  is  the  Wigner  of  is  the  Wigner  distribution  of  the  input 

u)  =  ^  J  f*(t-  \t)  f(t  +  It)  e~lTUJ  dr  (45) 

This  equation  can  be  seen  to  be  in  the  variables  t,uj,  and  is  in  general  a  partial  differential  equation 
of  twice  the  order  of  the  original  differential  equation,  that  is  order  2 n. 


4.1  Zero  Driving  Force  situation 

If  the  driving  force  is  zero  then 


P*(A)P(B)WX,X  =  0 


(46) 


17 


In  this  case  the  equation  of  evolution  can  be  simplified  to  two  equations 


P*(A)WX,X  =  P(B)WX,X  =  0  (47) 

These  equations  can  be  derived  in  two  equivalent  ways.  First,  directly  from  Eq.  (38)  with  /  =  0,  or 
one  can  start  from  the  definition  of  the  Wigner  distribution  and  show  it  straightforwardly.  We  note 
that  while  Eq.  (43)  is  a  differential  equation  of  order  2 n,  Eqs.  (47)  are  of  order  n. 


5  Ordinary  Differential  Equations  with  time  dependent  coefficients 
into  phase  space 

Consider  the  case  of  an  ordinary  differential  equation  with  time-varying  coefficients, 

,  .  dnx(t )  .dn~1x(t)  ,  dx(t ) 

Qn(*)  dfn  +On-i(t)  dtn_\  - \-ai(t)—^-  +  ao(t)x(t)  =  f(t)  (48) 

As  before  we  rewrite  this  in  polynomial  notation 

P(D,t)x(t)  =  f(t)  (49) 


where 

P{D,  t )  =  an(t)Dn  +  an-i{t)Dn~l . . .  +  a\(t)D  +  a0(t) 

A  similar  derivation  that  led  to  Eq.  (43)  leads  now  to 


P*(A1£)P(B,F)Wx,x(t,cu)  =  Wf,f(t,u) 

where 

„  1  d  —Id 

£  =  sa;  +  f  r=~nd^  +  t 

and  the  operators  A  and  B  are  given  in  Eq.  (44). 

For  the  zero  driving  force  case  one  obtains 


(50) 


(51) 


(52) 


P*(A,£)Wx,x(t,u;)  =  P(B,P)Wx,x(t,u ;)  =  0  (53) 

In  Eq.  (51)  P(B,P )  means  that  in  the  polynomial  P(D,t )  one  substitutes  B  for  D,and  T  for  t, 
and  P*(A,  S)  is  obtained  similarly.  Also,  the  complex  conjugation  in  P*(A,  £)  means  that  the  original 
polynomial  is  conjugated  and  not  the  arguments  A,  £  .  Note  that  Eq.  (51)  is  a  partial  differential 
equation,  and  that  makes  sense  since  we  are  dealing  with  two  variables  jointly,  namely  time  and 
frequency.  The  derivation  of  Eq.  (51)  is  given  in  the  Appendix. 
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6  Analysis  of  differential  Equations  in  Time- Frequency  Phase  Space 
for  the  Bilinear  Distributions 

We  consider  a  linear  system  defined  by  an  ordinary  differential  equation  of  the  type 

.  ,dnx(t)  ,dn~1x(t)  dx(t)  .  .  .  . 

On(t)  dfn  +  an-i(t)  dtn  \  ■■■  +  +  a0(t)z(t)  =  f(t)  (54) 

where  /(f)  is  the  input  or  forcing  term,  x(t)  is  the  output  of  the  system  and  an(t) , . . . ,  ao(t)  are 
the  time-varying  coefficients,  generally  complex.  As  mentioned  we  have  been  able  to  obtain  the 
governing  equation  for  the  Wigner  distribution  for  the  solution  of  this  equation.  We  consider  here 
the  possibility  of  transforming  such  equations  into  the  time-frequency  domain  for  a  general  bilinear 
class  of  distributions.  We  will  show  that  this  is  possible  for  a  bilinear  distribution  associated  to  the 
solution  of  Eq.  (54).  We  will  do  so  in  the  kernel  notation  and  the  general  K  formalism 


6.1  Equation  of  motion  in  the  kernel  notation 

The  general  class  of  bilinear  distributions  is 

C(t,u)  =  JJJ  x*(u  —  t/2)x(u  +  r/2)<j>c(9,  T)e~iet~lTUJ+ieudud9dT  (55) 

where  />c(0,t)  is  the  kernel.  To  obtain  an  equation  for  C(t,  u)  associated  to  Eq.  (54)  we  first  put  it 
into  polynomial  notation 


where,  as  usual 


and 


P(D,  t)x(t)  =  /(f) 


P{D ,  t)  —  an(t)Dn  +  an—i(t)Dtl  +  •  •  •  +  ai(t)D  +  ao(t) 


D  = 


dt 


(56) 

(57) 

(58) 


Then  the  equation  for  an  arbitrary  bilinear  time-frequency  distribution  C(t,  uj)  is 


P*(Ac,£c)P(Bc,Pc)Cx(t,u ;)  =  Cf(t,cj) 


(59) 


where 


Ac  = 

Bc  = 

£r  = 


1  d  1  d 
i  dt'  i  d lo 
1  8  1  d 
i  dt '  i  doj 
1  d  1  d 

i  dt '  i  dijj 

1  d  1  d 

i  dt '  i  duj 


1  d 


- ILU 


ILU 


2  dt 

1  d_ 

2  dt 


^  d_ 
^  2  idu 


t  — 


^d_ 

2  i  dco 


L~  1 


k- 1 


k  1 


1  d  1  d 
i  dt'  i  duj 
1  d  1  d 
i  dt'  i  du 
1  d  1  d 
i  dt'  i  du 
1  d  1  d 
i  dt'  i  du 


(60) 

(61) 

(62) 

(63) 
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We  have  already  derived  the  operators  Ac  and  Bc  in  [?,  ?].  Here  we  present  the  derivation  of  the 
operators  £c  and  Tc.  We  use  the  following  relation 

ca^)=-i>c(\§t,\l)wat,u)  (64) 

Now  we  evaluate 


°tf,g  =  <t>, 

=  <A 
=  <A 

=  £CCf,g 


1  d  1  d 

i  dt '  i  dtu 

~~)£W 
i  dt  i  olo 

1  d  1  d 
i  dt'  i  dio 


Wi 


tf,g 


f-g 


4-1 


1  d_ 

idt'idu)  f’9 


(65) 

(66) 

(67) 

(68) 


Where  we  have  used  that 


Wtf,g=£Wf,g 


Following  the  same  approach  one  has 


Cf,tg  ~  ^cCf^g 


We  can  now  evaluate 


(69) 

(70) 


(71) 


Ca{t)f,g  ~  CY,ant™f,g  (72) 

=  J2<Ctnf,g  (73) 

=  YJ<ZcCf,g  (74) 

=  a*(£c)Cftg  (75) 

With  the  identical  procedure  one  also  obtains 

Cm)g  =  b(Fc)Cf,g  (76) 


6.2  Equation  of  motion  in  the  K  notation 

All  the  bilinear  distributions  can  also  be  written  in  the  following  way 

C(t,  uj)  =  [[  Kc(t,u,t' x(t")  dt'dt" 


(77) 

j  J 

where  now  Kc(t,  uj,  t1 ,  t")  is  the  new  kernel  associated  to  every  distribution  C(t,u).  In  particular 
Kc(t,  uj,  t' ,  t")  is  given  by 

s(t,  uj,  x ’,  x)  =  J  x’  -  X)e-i{t-(x+x'Wed8  (78) 


Kr 
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We  now  claim  that  the  equation  for  the  K  notation  associated  to  Eq.  (54)  is 


P*(Ac,£c)P(Bc,Pc)Cx(t,u ;)  =  Cf{t,u>) 


(79) 


where  the  operators  have  the  same  form  given  in  Eqs.  (60)-(63),  where  one  has  to  make  the  substitution 


1  d  1  d 


1  d 


j.  is  j.  is  \  _  /  ~v~  I  .  ./  .  vy  ,  1  9  \  , , 

)  =  2?re  a“  /  ( *,&,*-*  +  -  7^—  ) e  at^ 


(80) 


2i(9tc  2  iduj/ 

To  prove  Eq.  (80)  one  considers  that  Eq.  (78)  is  basically  a  Fourier  transform  for  every  fixed  x,  x' . 
Hence  we  make  the  following  change  of  variables 


T  =  X  —  X 

t'  =  t  —  (x  +  x') /2 


(81) 

(82) 


which  gives 


x'  =  t  —  t'  T  t  2 

(83) 

x  =  t  —  t'  —  r/2 

(84) 

We  substitute  in 

Eq.  (78)  and  we  obtain 

J  (f>c(0,  r)e~lt  edd  =  4ir2etuJT  Kc(t,  u,t  —  t'  +  r/2,  t  —  t'  —  r/2) 

(85) 

Then  we  obtain 

<t>c(9,  t )  =  2irelu>T  J  Kc(t,  u,t  —  t'  +  r/2,  t  —  t'  —  r /2)elt  9 dt' 

(86) 

By  substituting 

1  3 

(87) 

1  <9 

i 

(88) 

we  finally  obtain  Eq.  (80). 


7  Differential  equations  for  the  Short-Time  Fourier  Transform  do¬ 


main 


We  define  the  Short-Time  Fourier  transform,  Sx(t,tv),  of  a  signal,  x(t),  by 

Sx(t,uj)  =  _  f  h(r  —  t)x{T)e~lru) dr 

v27t  J 


(89) 


where  h(t)  is  the  window  function.  The  equation  for  the  Short-Time  Fourier  Transform  associated  to 
Eq.  (54)  is 

P{As,£s)Sx(t,<jj)  =  Sf(t,u)  (90) 
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where 


As  —  \  —  +  iw 


£s  =  i~ 


d_ 
Of, 
,0_ 
Ou 


To  prove  Eq.  (90)  we  first  evaluate 


u)  =  --L  f  h{T-t)^^e~iTUJdT 

y/2 n  J  dr 


1 


\/2tt 


(h(T-t)x(T)e-^)+_Z  + 


1  f  f  dh[r  —  t ) 


V2^ 


dr 


x(r)e  lTU>  —  iu)h{r  —  t)x(r)e  ZTU>  ]  dr 


^  h(r  —  t)x(r)e  lTiJ dr  +  iuj  _  f  h(r  —  t)x(r)e  lTUdr 


\Z2tt  Of 


y/2 7T  . 


=  ^  +  iu^j  Sx{t,LO ) 

=  Assx(t,uj ) 


where  we  have  used 


(91) 

(92) 

(93) 

(94) 

(95) 

(96) 

(97) 

(98) 


Then  we  consider 


h(— oo)  =  h(+ oo)  =  0 
dh(r  —  t)  Oh(r  —  t ) 

dr  Of 


j  w) 


h(r  —  t)a(r)x(r)e  *ra,<iT 


=  a{it)7^Ih{T-tMr)e"T“dT 

=  a(£s)Sx(t,u ;) 


(99) 

(100) 

(101) 

(102) 

(103) 


which  gives  Eq.  (90). 


8  Transformation  to  the  Wavelet  domain 


We  now  consider  the  problem  of  writing  an  equation  for  the  Continuous  Wavelet  Transform  (CWT) 
whose  solution  corresponds  to  the  time  domain  function  x(t)  as  given  by  Eq.  (54).  Here  we  present  a 
restricted  version  of  the  problem,  addressing  linear  differential  equations  with  constant  coefficients 


dnx(t )  dn  1x(t) 

,n^^  +  Cln-1  dt "-i 


•  •  •  +  a\ 


dx(t) 

dt 


+  a0x(t) 


m 


(104) 


22 


The  CWT  of  a  signal  x(t)  is  defined  as 

Cx(a,b )  =  ~^=  J  ip*  x{t)dt 

where  a  >  0,  — oo  <  b  <  +oo.  We  rewrite  Eq.  (104)  in  the  simplified  polynomial  notation 

p(D)X(t)  =  m 

The  equation  for  the  CWT  is 


where 


To  prove  Eq.  (107)  we  evaluate 


P(Aw)Cx(a,b)  =  Cf(a,b) 


A  -1 

w  db 


1 

y/a 

/’ 

1  f 

yfd 

1  V 

1 

d 

t  —  b\  dx(t) 


a  J  dt 


dt 


t  —  b 


x{t) 


+00 


j=  (— i  ^ 

\/a  ./  at  \  a 


\fa  db  J  1 
=  AwCx(a,,  b) 


t  —  b\  dx(t) 


a  J  dt 


dt 


where  we  have  used 


'ip(-oo)  =  V^+oo)  =  0 
d  ft-b\  d  ft-b 

It*  —  r-m*  (  — 


(105) 

(106) 

(107) 

(108) 

(109) 

(110) 

(111) 

(112) 

(113) 

(114) 


Also,  we  note  that  since  Cx(a,b)  is  a  linear  operation,  then,  for  any  given  complex  constant,  k,  we 
have  that  Chx{a,  b)  =  kCx(a,  b).  This  fact  was  used  in  the  above  derivation 


9  Transformation  to  the  Ambiguity  function  domain 

We  have  considered  a  number  of  topics  that  involve  two  fundamental  problems  in  signal  analysis,  sys¬ 
tems  governed  by  differential  equations  and  whose  solutions  generally  give  a  nonstationary  spectrum. 
We  also  point  out  that  while  in  this  paper  we  have  just  considered  ordinary  differential  equations 
these  methods  can  be  readily  extended  to  situations  governed  by  partial  differential  equations,  and 
to  systems  of  ordinary  differential  equations.  In  conclusion  we  also  mention  that  we  have  worked 
out  differential  equations  for  the  ambiguity  function  and  the  autocorrelation  function.  The  ambiguity 
function  is  defined  by 

Ax(t,u)  =  —  [  x*{t  —  r/2)x(t  +  T/2)el6tdt 

27 r  / 
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(115) 


and  the  differential  equations  that  is  satisfied  is 


P*(Aa,£a)P(Ba,Jra)Ax(t,U})  =  Af(t,Cd) 


where 


2  dr 

id  1 
2r’ 


For  the  autocorrelation  function,  defined  by 


Ba-~\i9+^ 

_  id  1 

Pa  —  ~  WX  +  XT 

i  d0  2 


Rx{t1,t2)  =  E[X{t1)X\t2)\ 


the  differential  equation  is  given  by 


(116) 

(117) 

(118) 

(119) 

(120) 


10  Partial  Differential  Equations 


For  partial  differential  equations  of  the  form 

N  M 

'  dxk 


dk  1V1  dl 

^2ak(x,t)^u(x,t) +  ^2bi(x,t)-^u(x,t)  =  f(x,t ) 


(121) 


k= 1  1=1 

it  is  possible  to  obtain  an  equation  in  the  Wigner  distribution  domain  similar  to  the  ordinary  differential 
equation  case,  provided  that  we  use  the  four  dimensional  Wigner  distribution,  Z.  defined  as 

1 


Zu,u(x i  Pjt,  Ld)  — 


(20 


u*(x  —  \rx ,  t  —  \t)u(x  +  \tx,  t  +  \r)e 


1  „ — irco — 


lTxPdrdTr 


(122) 


This  distribution  is  discussed  in  the  Appendix,  here  we  have  reported  its  definition  again  for  clarity. 
We  shall  need  to  use  the  following  operators, 


1  d_ 

2  dx 
1  d 


Ax  =  -  — - in  At  =  -— - lid 

°  Q--  y  2  dt 


1  d_ 

2  dt 
i  d 


^=0^:  +  ^  Bt=  2  dt  +iuJ 


and 


2  dx 
1  d 


1  d 


£-X  —  „  .  o  P  X  J~X  —  P  x 


£t  = 


2  j  dp 
^  d_ 

2  j  did 


2  j  dp 
_  Id 

+  t  Bt  =  ~Yjdu+t 


(123) 

(124) 

(125) 

(126) 


The  same  approach  of  the  previous  sections  for  ordinary  differential  equations  gives 


N 


M 


Y/a,i(£x,£t)Ak  +  Y/b*l(£x,£t)A 


Lfc=l 


1=1 


N 


M 


Y  am(Bx,Bt)B™  +  Y  bn(Px,Pt)B, 


m=  1 


n=  1 


Zu,u=  ZfJ  (127) 
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If  the  driving  force  f(x,t)  is  zero  one  has  that 


N  M 

E  [a%(ex,£t)Ai\  Zu,u  =  -E  [bm,St)A 


k=  1 
N 


E  [*k{Fx,Ft)Bkx 


1=1 

M 


Zn, 


ZUltt  =  -E 

fc=i  “  “  i=i 

We  have  found  it  convenient  to  combine  these  two  equations  by  adding  and  subtracting  them 


(128) 

(129) 


N  M 

E  [a%{Sx,St)Akx  ±  ak(^x,^t)Bk]  ZU)U  =  E  [&?(£*,  £t)A\  ±  6,(^X)  (130) 


fc=i 


/=! 


10.1  Notation  and  Definitions 

All  the  integrals  without  limits  mean  integration  from  — oo  to  oo. 

Ordinary  differentiation  of  functions  with  respect  to  time  will  be  indicated  by 

fj  fjn 

m  =  Jtg(t)  9<"'  =  (i3i) 

We  define  the  Wigner  distribution  for  a  signal  x{t)  by 

WXjX{t,uj)  = -^  j  x*  (t  -  \t)  x(t  +  \t)  e~lTUJ  dr  (132) 

and  the  cross- Wigner  distribution  between  two  signals,  x(t)  and  y(t)  by 

WXuX2(t,u)  =  ^-  f  x\{t-  \t)  x2(t  +  \t)  e~lTU  dr  (133) 

When  we  deal  with  partial  differential  equations,  we  need  to  consider  multidimensional  signals 
(fields),  u(x,t).  For  a  field  the  Wigner  distribution  is 

WUtU(x,p,t )  =  ^  J  u*(x-  \rx,t)u{x+  \Tx,t)e~lTxV  drx  (134) 

and  analogously  the  cross-Wigner  distribution  between  two  fields,  u\{x,t)  and  u2(x,t)  is 

Wul,U2{x,P,t)  =  J  u\(x  -  \rx,t)u2(x  +  \Tx,t)e~lTxP  drx  (135) 
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For  partial  differential  equations  it  is  generally  not  possible  to  write  an  equation  for  the  Wigner 
distribution,  WUiU(x,p,  t),  corresponding  to  an  arbitrary  equation  governing  the  field.  It  is  nevertheless 
possible  to  always  derive  such  an  equation  for  the  more  general  Wigner  distribution  which  we  define 
by 


Zu,u{x,P,t,uj)  =  *  2  J  u*(x-\Tx,t-\T)u(x+\Tx,t+\T)e  tTUJ  lTxP  dr  drx  (136) 

We  note  that  the  ordinary  Wigner  distribution,  W(x,p,t),  may  be  obtained  from  Z(x,p,t,u>)  by  way 
of 

WUjU(x,p,t)  =  J  ZUtU(x,p,t,uj)duj  (137) 

A  significant  simplification  in  notation  is  achieved  by  defining  the  following  operators 


.  Id. 

A‘=2  9t~'u 

_  Id 

Bi  =  2  U+,U 

(138) 

c  1  9 

t  ^  2i  dco 

Si 

II 

CH* 

1 

(139) 

Operators  of  this  kind  were 

defined  by  Moyal.  When  dealing  with  a  field  u(x, 

t)  one  has  to  introduce 

the  additional  operators 

a  id. 

Ax  =  -  — - ip 

2  dx  F 

D  19 

b-  =  2  ai  +  ,p 

(140) 

„  1  d 

£*  =  X+2i¥p 

r  i  d 

Tx  X  2idp 

(141) 

Also,  we  will  indicate  ordinary  differentiation  in  the  following  alternative  ways 

d  dn 

m  =  Jtg(t)  9«">  =  ^s(t)  (142) 

and  we  will  use  the  differential  operator, 

D  =  1  <143) 

As  mentioned  above,  we  will  see  that  it  is  not  always  possible  to  obtain  an  equation  of  motion  for 
the  Wigner  distribution  WUjU(x,p,t),  but  it  is  always  possible  to  obtain  an  equation  of  motion  for  a 
Wigner  type  distribution  of  the  four  variables,  namely  position,  momentum  time  and  frequency.  We 
define  such  a  function  by 

Zu,u(x,P,t,uj)  =  2  J  u*  (x  -  \rx,t  -  \t)u{x  +  \rXlt  +  \T)e~lTU1~lTxP  dr  drx  (144) 
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The  relation  between  W  and  Z  is  that  W  is  the  marginal  of  Z 


J  ZUjU(x,p,  t,  u)  dio  =  WUtU(x,p,  t) 


It  is  convenient  to  define  the  following  operators1 

.  Id. 

=  -  —  -tu 


2  m 


< St  —  t-\- 


^ a 

2  i  duj 


D  19 

Bt  =  V  +  IV 


2  dt 


Ft  =  t  — 


1_ J) 

2  i  duj 


(145) 


(146) 

(147) 


Also,  the  operators  Ax,  BXl  Ex,  Fx  for  fields  are  similarly  defined  as  per  Eq.  (146)  and  Eq.  (147), 
where  t  is  replaced  by  x  and  lu  by  the  momentum,  p. 


11  Example:  x  =  f{t) 


We  consider  an  example  that  to  some  extent  is  analytically  doable.  For  the  differential  equation  we 
take 


x  =  fit ) 


with 


fit)  =  cie-^2/2+i/3*72 

which  with  the  boundary  condition  that  .x(— oo)  =  0  gives 


=  C  c-at2 /2+iflt2 /2 

a  —  i/3 


x(t)  = - —e 


(148) 

(149) 

(150) 


The  Wigner  distribution  of  both  x  and  /  can  be  done  exactly 

u)  =  2A  fZ  [  at2  +  {u_  m)2/a  _  .  ] 


1  I  IT  C2 


0—at2  —  (u)—/3t)2 /a 


(151) 


7T  V  a  a2  +  /32 

We  now  show  how  our  procedure  can  be  used  to  obtain  the  Wigner  distribution  directly.  Using 
Eq.  (39)  the  equation  of  motion  for  the  Wigner  distribution  is 


that  is 


ABWXtX{t,u)  =  Wfj(t,v) 


1  d 

Vx,xit,u)  +  u2Wx,x(t,u>)  =  Wfjit,uj ) 


(152) 


(153) 


Operators  of  this  kind  were  defined  by  Moyal. 
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This  equation  is  a  partial  differential  equation  but  since  there  are  no  derivatives  with  respect  to 
u  we  can  solve  it  as  an  ordinary  differential  equation  keeping  to  fixed.  By  repeating  the  process  for 
a  range  of  oj  we  hence  get  a  numerical  solution  for  Wx  x(t,(jS).  We  have  done  so  and  compared  to 
the  exact  answer,  Eq.  (151),  with  the  result  that  the  answers  are  identical  to  within  a  fraction  of  a 
percent. 


12  Example:  RC  circuit  with  deterministic  chirp  driving  force 


As  an  example  consider  an  RC  circuit  that  is  being  driven  by  a  driving  force  whose  instantaneous 
frequency  varies  linearly  with  time.  The  governing  equation  is 

d^l  +  kx  =  eiU0t+i^2l2  (154) 

at 

where  A;  is  a  positive  constant.  Using  Eq.  (43)  we  have  that  the  associated  equation  for  the  Wigner 
distribution  is 

[A  +  k\[B  +  k]WXtX(t,u})  =  5(uj  —  loq  —  fit)  (155) 


This  equation  can  be  solved  exactly,  giving 


e-2kr 


sin  2c or 
2  ui 


where  u(t)  is  the  step  function  and 


As  initial  conditions  we  have  taken 


T  =  t  — 


CO  —  CUo 

“ d~ 


(156) 


(157) 


WX:  x(-oo,u)  =  0 
d 

-q-Wx,x{-oo,uS)  =  0 

that  can  be  proved  to  correspond  to  the  case 


(158) 

(159) 


x(— oo)  =  0 

In  Fig.  13  we  give  a  specific  case  where  k  =  20,  ujq  =  0,  and  (3  =  113.  The  dashed  line  represents  the 
instantaneous  frequency  of  the  input  chirp  f(t), 

uji  =  +  /3t  =  (3t  (160) 

while  the  gray  image  is  the  Wigner  distribution  obtained  from  Eq.  (156)  with  the  given  parameters. 
The  figure  clearly  reveals  the  bandpass  behavior  of  the  system,  that  gradually  filters  out  the  forcing 
term  /(f). 


28 


t(s) 


Figure  13:  The  dashed  line  indicates  the  instantaneous  frequency  of  the  input  chirp.  The  gray  image 
is  the  Wigner  distribution  WXjX(t,uj)  given  in  Eq.  (156)  and  corresponding  to  the  solution  x(t)  to 
Eq.  (154).  The  Wigner  distribution  highlights  the  lowpass  behavior  of  the  system,  that  filters  out  the 
input  signal  as  t  — >  oo. 

13  Example:  RLC  circuit  with  generic  deterministic  input 

Consider  the  harmonic  oscillator  with  a  deterministic  driving  force, 

^+^+ulx-m  (!«) 

where  x(t)  is  the  state  variable  (e.g.  current,  position)  and  f(t)  is  the  driving  term.  If  we  want  to 
study  the  time-frequency  properties  we  could  solve  this  equation  and  substitute  the  answer  into  the 
Wigner  distribution, 


Wi,i(i,w)  =  ^  /  x*(t  -  \t)  x(t  +  \t)  e  lTUJ  dr  (162) 

We  want  to  write  the  equation  of  motion  for  the  Wigner  distribution  and  solve  that  directly.  We  first 
rewrite  the  equation  in  polynomial  notation 

[  D2 +  2^D  +  uj'q]  x(t)  =  f(t)  (163) 

Using  our  method  for  ordinary  differential  equations  with  constant  coefficients,  Eq.  (43),  we  obtain 
the  equation  of  motion  for  the  Wigner  distribution 

[  +  2 fiA  +  Uq  ]  [  B2  +  2/iE>  +  lJq  ]  Wx,x  =  Wftf  (164) 

Explicitly, 

-  q4  q3  q2  q 

a*dB+a3W  +  a2w+aim+a° 


WXjX(t,u>)  =  io)  (165) 
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where. 


a  o  =  (cug  —  cp2)2  +  Afi2u j2 

(166) 

CL\  —  2/t((Uq  -f"  (U2) 

(167) 

°2  =  2^0  +  ^  +  2/i2) 

(168) 

G3=2M 

(169) 

04  =  1/16 

(170) 

14  Example:  The  Exact  Solution  to  the  Gliding  Tone  Problem 

In  1948  Barber  and  Ursell  [1]  and  independently  Hok  [3]  considered  the  problem  of  the  response  of 
a  harmonic  oscillator  to  a  “gliding  tone”.2  Specifically  the  issue  is  the  behavior  of  the  solution  to  a 
resonant  circuit  [1,3,4] 

+  2 74^  +  u2x(t)  =  eMt+W  (171) 

with 

f(t)  =  ei/3t^2  (172) 

Subsequent  to  Barber  and  Ursell,  and  Hok,  many  investigators  have  considered  this  problem  in  a 
variety  of  contexts  and  have  tried  to  qualitatively  understand  the  solution  and  also  obtain  approximate 
solutions.  An  exact  solution  to  this  problem  has  not  been  achieved.  We,  also,  have  not  been  able  to 
obtain  an  exact  explicit  solution;  but  we  have  been  able  to  obtain  the  exact  solution  to  the  Wigner 
distribution  of  x(t)l  We  have  been  able  to  obtain  the  exact  solution  by  using  our  method,  that  is  by 
transforming  the  original  equation  in  time  to  the  domain  of  the  Wigner  distribution,  and  by  solving 
it. 

We  give  here  the  explicit  solution  to  the  gliding  tone  problem  and  subsequently  we  give  a  few 
numerical  examples.  The  reason  this  is  called  the  gliding  tone  problem  is  because  the  instantaneous 
frequency  of  the  driving  force  increases  linearly, 

=  loi  +  (3t  (173) 

In  the  gliding  tone  problem  one  wants  to  ascertain  the  instantaneous  frequency  of  the  response.  There 
have  been  a  number  of  studies  made  by  examining  approximate  solutions  of  Eq.  (171),  because  indeed 
an  exact  solution  to  Eq.  (161)  with  /(t)  given  by  Eq.  (171)  has  not  been  achieved.  However,  we  have 
been  able  to  solve  the  equation  for  the  Wigner  distribution  of  the  gliding  tone  problem  exactly.  That 
2The  phrase  “gliding  tone”  was  used  by  Barber  and  Ursell. 
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Figure  14:  Wigner  distribution  of  the  solution  to  the  gliding  tone  problem.  Underdamped  case  with 
H  =  0.5. 

is,  we  have  been  able  to  solve  Eq.  (165)  when  the  input  term  is  the  Wigner  distribution  of  the  gliding 
tone,  Eq.  (172).  The  detailed  answer  is  given  in  Sect.  14.1. 

We  now  give  some  graphical  examples  to  illustrate  the  results.  We  first  consider  the  underdamped 
case,  that  is,  when  //  <  loq.  In  Figs.  14-16  we  plot  the  Wigner  WXjX(t,uj)  for  the  three  cases  /i  = 
0.5,  fi  =  1,  and  /i  =  1.5,  respectively  .  In  both  of  the  three  cases  we  chose  luq  =  18  and  uj\  =  (3  =  4.  The 
gray  scale  image  in  every  picture  is  the  exact  Wigner  distribution  of  x(t);  the  dashed  line  represents 
the  instantaneous  frequency  of  the  forcing  chirp,  that  is  uq(f)  =  uq  +  f3t.  The  input  chirp  /(f)  is 
concentrated  only  along  this  line,  because  its  representation  in  the  Wigner  distribution  domain  is 
S(u  —  uj\  —  /3t).  We  see  that  the  response  of  the  system  is  mainly  concentrated  around  the  critical 
frequency  toc  (it  is  u)c  ~  wo)>  while  it  is  weaker  at  all  the  other  frequencies.  Also,  observing  the  limit 
at  lo  =  luc,  one  can  see  that  the  Wigner  distribution  has  an  exponential  damping  factor,  where  the 
damping  coefficient  is  2/i,  which  is  twice  the  damping  of  the  free  oscillation  factor  [i  of  the  system. 
Comparing  the  three  pictures  we  see  how  changing  the  damping  factor  [i  influences  the  system  response. 
Smaller  values  of  q  imply  less  damping  and  hence  longer  tails  in  the  response  along  loc.  Increasing 
/i  forces  the  system  to  have  a  stronger  damping  and  that  is  reflected  in  the  shorter  tail  of  the  main 
response  located  around  the  resonant  frequency  ojc. 

In  Fig.  17  we  give  an  example  of  an  overdamped  case  where  q  >  loq,  and  in  particular  we  take 
H  =  30.  Here  the  system  response  is  Inharmonic,  and  we  do  not  have  any  special  resonant  frequency. 
Notice  that  the  output  is  greater  for  small  times  f,  while  when  t  — >  oo  the  response  goes  to  zero.  This 
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Figure  15:  Wigner  distribution  of  the  solution  to  the  gliding  tone  problem.  Underdamped  case  with 
I1  =  1 


Figure  16:  Wigner  distribution  of  the  solution  to  the  gliding  tone  problem.  Underdamped  case  with 
//  =  1.5 
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Figure  17:  Wigner  distribution  of  the  solution  to  the  gliding  tone  problem  for  a  overdamped  case, 
M  >  ^o- 


Figure  18:  Wigner  distribution  of  the  solution  to  the  gliding  tone  problem  for  a  critically  damped 
case,  fjj  —  cjq* 
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is  in  complete  agreement  with  the  result  obtained  considering  the  system  transfer  function. 

Finally  in  Fig.  18  we  show  a  critically  damped  case,  with  //  =  cuc  =  18.  Considerations  on  this 
case  are  similar  to  those  for  the  overdamped  case. 

We  point  out  that  we  have  obtained  the  Wigner  distribution  from  Eq.  (165)  by  choosing  the 
following  initial  conditions: 


dWx,x(- oo,cu)  d2Wx,x(-oo,u)  d3Wx,x(-oo,u) 

OO^iO)  — 


dt 


dt2 


dt3 


=  0 


(174) 


We  have  proven  XciteRepresent  ability  that  this  choice  corresponds  to  finding  the  Wigner  distri¬ 
bution  of  the  solution  x(t )  that  has  zero  initial  conditions  at  t  =  — oo. 


14.1  The  exact  solution  to  the  Gliding  Tone  Problem 

We  now  give  the  exact  solution  of  the  Wigner  distribution  of  x{t)  which  satisfies  Eq.  (165) 


Z2  ~  Z 1 


Zl  -  Z 1 


e-2zi r  _  e-2z2 r  g-2zir  —  e~2z2T 


Z2  ~  Z 1 


Z2  ~  Z 1 


+ 


g-2z2T  _  g— 2z2"t  g-2zir  _  g-2z2T 


Zl  ~  Z2  \  Z 2  -  Z2 


Z2  ~  Z 1 


with 


and 


T  =  t  —  Uj/ /3 


Zi  =  -ju  +  n  -  \/n 2  -  wg 


zi=ju  +  /I  -  \/  n2  -wg 


z2  =  -juj  +  n  +  \  nz 


Z2=JU  +  li  +  -  W{ 

and  where  u(t)  is  the  Heaviside  step  function  given  by 

1,  t  >  0 


(175) 

(176) 

(177) 

(178) 

(179) 

(180) 


u(t)  =  { 


(181) 


0,  t<  0 


14.1.1  Underdamped,  Overdamped,  and  Critically  Damped  Cases 

We  now  explicitly  specialize  to  the  underdamped,  overdamped  and  critically  damped  cases.  As  is 
standard  we  define  the  critical  frequency,  loc,  for  these  three  cases 
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uic=  \Juj%-  p2 

P<U  o 

Underdamped 

Uc=  \J p2  -u>l 

/i  >  cu0 

Overdamped 

wc  =  0 

n  =  u  0 

Critically  damped 

The  explicit  Wigner  distributions  are 
Underdamped: 


W{t^)  =  %bcT)e^x 


sin(2(w  —  ujc)t)  sin(2(w  +  uic)t) 


Overdamped: 


W(t,u)  =  y-jy u(r)e  2/ir  x 


w(<W  —  Wc)  L)(U)  +  Wc) 

sin(2wr)  cosh(2wcr)  cos(2wr)  sinh(2wcr) 

w(w2+w2)  wc(w2  +  w2) 


Critically  Damped: 


WM  =  ±-u(T)e-^Sin{2lJT)  -2J C°S(2"T| 


(182) 


In  the  above  solutions  there  are  singularities  at  some  values  of  w.  We  give  the  limits  at  those 
singular  values  for  the  three  cases: 

Underdamped: 

1  ,  ,  [4u>r.T  —  sin (Aur.r)  1  ,  . 

(Too) 


lim  W{t,u)  = ——u{T)e  2^r 

LO^±UJc  2  \jj\UJc 


4  loct  —  sin(4wcr) 


2cu2 


Overdamped  : 


Critically  Damped: 


lim  W(t,  u)  =  |^T u{t)e  2/it 


lim  W(t,ui)  =  ^—u(t)e  2 pt  x 

\(j\ 


sin(2wct)  —  2wcf  cos(2rocf) 


OJt 


2uct  cosh(2wct)  —  sinh(2wcf) 


OJ'i 


(184) 


(185) 


15  Example:  a|^  +  V(x,t)u(x,t)  = 

To  bring  forth  some  of  the  ideas  discussed  above,  let  us  consider  the  following  type  of  equation 

c)^ii  chi 

a~7^2  +V(x,t)u(x,t)  =  b—  (186) 
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where  for  the  moment  we  do  not  make  any  assumptions  about  a  or  b,  except  that  they  are  com¬ 
plex  constants.  This  leads  to  the  following  equations  for  the  Wigner  distribution  for  fields  (see  the 
Appendix) 


a*A2xZu,u  +  V*(£x,£t)Zu,u  =  b*AtZu, v 
aB2xZu,u  +  x,Ft)Zu,u  =  bBtZU)U 
Adding  and  subtracting  the  above  two  equations  we  have 


(187) 

(188) 


[a*A2x  ±  aBl ]  ZUyU  +  [V*(£x,£t)  ±  V{Fx,Ft))  Zu,u  =  [b*At  ±  bBt ]  Zu,u 

The  case  where  b  is  pure  imaginary  and  a  real 

We  consider  the  special  case  where  a  is  real  and  b  pure  imaginary  and  write 

b  =  i\b\ 

Taking  the  negative  sign  in  Eq.  (189),  we  have 

a  [Al  -  Bl]  Zu,u  +  [V*(£x,£t)  -  V(Tx,Tt)]  Zu,u  =  -i\b\  [At  +  Bt]  ZU>1 


But 


Al~Bl  =  -2ip 


At  +  Bt  = 


d_ 

dt 


and  therefore 


or 


dZ  d 

-i\b\-^£  =  -i2ap—ZU)U  +  [V*(£x,£t)  -  V(Tx,Tt)}  Zu,v 


3Z  d 

\b\  =  2 ap—Zu>u  +  i[V*{£x,£t)  -  V(Tx,Tt)\Zu,v 


(189) 


(190) 

(191) 

(192) 

(193) 

(194) 

(195) 


This  is  the  equation  of  motion  in  its  most  general  terms.  It  can  be  written  in  a  more  explicit  way. 
First  we  mention  that  for  any  four  dimensional  function  K(x,p,t,  to) 


V(J7x,J7t)K(x,p,t,u)  = 


Z'j  j  V(x  +  K(x,p' ,t,u')dx'dp'dt'dw'  (196) 


and 


V*(£x,  £t)K(x,p,  t,  u)  = 


1 

2vr 


J  V*(x  +  x'/2,t  +  K(x,p\t,J)dx' dp' dt'dJ  (197) 
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Now  consider  real  potentials,  then 


[V(£x,£t)  -  V(Fx,Ft)\  ZUiU(x,p,t,u)  =  J  V(x  +  x'/2,t  +  t'/2)x 

sin[(//  —  p)x'  +  ( u '  —  u)t']ZUtU(x,p' ,  t,  u')dx'dp'dt'du' 

and  hence  we  have  for  the  equation  of  motion 

r)7  r)  2  f 

\b\ =  2ap-^ZUiU  +  sin[(p'  -  p)x'  +  (u/  -  c o)t'}  x 

V(x  +  x' /2,  t  +  t' /2)ZUjU(x,p' ,  t,  u')dx' dp' dt' du'  (198) 

which  can  also  be  written  as 

c)7  f)  q  r 

\b\ =  2ap—Zu,u  +  ^2  /  sin[2(p'  -  p){x'  -  x)  +  2{u'  -  u)(t'  -  f)]  x 

V ( x t')ZUtU(x,  p  ,  t,  u')dx'  dp'  dt'  du'  (199) 

15.1  The  Schroeder  Equation 

Let  us  now  specialize  explicitly  to  the  Schroeder  equation. 

|6|  =  1  °=~  (200) 

By  using  the  results  above  we  have 

^ft4  =  -mFq2*’*  +  7  [V‘i£“’£,)  "  Z**  <201) 

and 

dZ^___p_d_  8 

dt  mdq  ^  (27t)2 

J  sin[2(j/  —  p)(q  —  q)  +  2(u'  —  u)(t'  —  t)]V(q  ,t')Z^^(q,p'  ,t,u')dq'dp'dt'du'  (202) 

These  are  the  equations  of  motion  for  the  Wigner  distribution  for  time  dependent  potentials.  Notice 
that  one  cannot  integrate  out  u  and  get  a  differential  equation  for  the  ordinary  Wigner  distribution. 
However,  if  we  have  time-independent  potentials  then  one  can  integrate  out  u  to  obtain  such  an 
equation.  We  take 

V(q,t)  =  V(q)  (203) 

and  use 

J  sin[(p7  —  p)q  +  ( u '  —  u)t']dt'  =  2tt  5{u  —  u')  sin[(p7  —  p)q]  (204) 

then  Eq.  (202)  reduces  to 
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dWvv0  =  +  1  f  V(q  + q'/2)sin[(p' -  p)q/]Wi)^(q,p',t)dq/dp'  (205) 

or 

=  ~m^qW^  +  nj  sin[2^  ~  ~  q)]w^A^P^  t)d<fdp'  (206) 

which  is  the  equation  originally  obtained  by  Wigner  and  Moyal. 

Liouville  sine  operator  form  of  the  equation  of  motion.  Moyal  derived  the  equation  of  motion 
for  the  Wigner  distribution  in  a  Liouville  type  operator  form.  This  was  done  for  time-independent 
potentials.  We  now  ask  whether  the  same  type  of  equation  can  be  obtained  for  the  time-dependent 
case.  First  we  recollect  the  Liouville  type  expression  Moyal  XciteMoyal  obtained  for  Eq.  (206).  We 
shall  present  the  result  in  a  somewhat  different  way  then  Moyal.  The  main  formula  to  be  used  is  that 
for  any  function  of  q,K(q),  and  any  function  of  q.p  .  M(q,p)  one  has  that  3 


1 

sin  - 
2 


d  d  d  d 


dpK  dqM  dpM  dqK 


K{q)M(q,p)  = 


TT 


sin[2(p'  —  p)(q'  —  q)]K(q,)M(q,p,)dp'dq'  (207) 


Therefore  Eq.  (206)  can  be  put  in  the  form 


dW^  p  d  .  1 

= - w-W^^,  +  2sin- 

moq  2 


dt 


d  d  8  d 


dpw  dqy  dpv  dqw 


V(q)WM(q,p,t) 


(208) 


But  we  also  have  (trivially)  that 


p  d  TTr  l 

- =  2sm  n 

moq  2 

**>  =  k 


d  d  d  d 


dpw  dqA  dp  a  dqw 


A(p)W^(q,p,  t ) 


(209) 

(210) 


and  therefore  if  one  defines 


H(q'p) = L + v(q) 


then  we  have  the  Moyal  form  of  the  equation  of  motion, 


(211) 


dW^Tf,  .  1 

— =  2  sm  - 
dt  2 


d  d 


d  d 


dpw  dqH  dpu  dqw 


H(q,p)W^(q,p,  t ) 


(212) 


3 To  change  the  sign  of  the  right  member  to  a  plus,  just  switch  the  order  of  the  differentials  in  the  sine  bracket. 
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16  Example  §  = 


As  with  the  case  of  ordinary  differential  equations  we  will  illustrate  our  method  with  specific  examples. 
The  examples  we  use  are  the  Schroeder  free  particle  equation  and  the  diffusion  equation 


difi  .  d 

dt  a  dx 2 

dii  _  nd2u 
dt  dx2 


Schroedinger  free  particle;  a  =  h/(2m) 
Diffusion  equation;  D  =  diffusion  coefficient 


(213) 

(214) 


We  have  chosen  these  equations  because  they  are  fundamental.  In  addition  they  are  superficially 
similar.  We  want  to  illustrate  how  these  equations  compare  in  the  Wigner  representation,  that  is,  in 
phase  space.  The  Wigner  distribution  for  the  Schroeder  equation  has  been  studied  for  over  70  years 
but  we  believe  it  has  not  been  applied  to  the  diffusion  equation  and  equations  of  that  type.  The 
respective  momentum  functions  defined  by 


4>{p,  t) 
U(p,t ) 

satisfy  the  following  equations  of  motion 


i/j(x,t)e  lxpdx 
u(x,  t)e~lxpdx 


dip 

dt 

dU 

~dt 


=  — iap2(j) 
=  —Dp2U 


(215) 

(216) 


(217) 

(218) 


The  Wigner  distribution  combines  both  representations,  that  is,  it  is  a  joint  representation  of  position 
and  momentum.  However,  we  point  out  a  fundamental  difference  between  the  interpretation  of  the 
solution  of  the  Schroeder  and  diffusion  equation.  In  the  case  of  the  Schroeder  equation,  \ip(x,  t)\2 
and  \cj>(p,  t)\2  are  the  densities  of  position  and  momentum.  However,  in  the  case  of  diffusion,  u(x,t) 
and  U(p,t)  are  the  densities.  Thus,  in  the  case  of  the  Schroeder  equation  the  Wigner  distribution 
satisfies  the  so  called  marginal  conditions,  but  that  is  not  the  case  for  the  diffusion  equation.  None  the 
less  the  Wigner  distribution  gives  an  indication  on  how  momentum  and  position  are  jointly  related. 
More  precisely  one  should  think  of  the  representation  as  a  joint  representation  in  position  and  spatial 
frequency. 

We  recall  the  Wigner  distribution  for  a  field,  u(x,  t ) 

Wu(x,p,  t)  =  J  u*(x  -  ^t,  t)  u(x  +  ^t,  t)e~lTpdr  (219) 

=  ^J  U*(p+h,t)U(p--2e,t)e-iexde  (220) 

We  will  use  i/j  and  u  to  signify  the  solution  of  Schroeder  and  diffusion  equations  respectively  and  use 
W.xp(x,p,t )  and  Wu(x, p,t)  for  their  respective  Wigner  distributions.  When  we  apply  our  method  to 
the  above  equations  we  obtain 
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(221) 


dW,i, 

dt 

dWu 

dt 


—2 pa 


dWp 

dx 


Dd2Wu 
2  dx2 


2Dp2Wu 


(222) 


Eq.  (221)  was  first  obtained  by  Wigner  and  Moyal  and  its  properties  have  been  studied  for  many 
years. 

It  is  quite  interesting  that,  while  the  only  difference  between  the  two  original  equations,  Eq.  (213) 
and  Eq.  (214),  is  an  i,  the  difference  in  the  Wigner  distribution  equation  of  motion  is  quite  dramatic. 
More  importantly  we  will  see  that  in  the  Wigner  domain  both  the  mathematics  and  insight  become 
clearer.  We  point  out  that  the  results  we  present  for  the  Schroeder  equation  are  classic  in  the  work  of 
Wigner,  Moyal,  and  many  others,  but  the  results  we  present  for  the  diffusion  equation  we  believe  to 
be  new.  We  mention  that  these  equations  may  be  related  to  the  respective  Fokker-Planck  equations 
but  that  will  not  be  pursued  here. 

We  point  out  that  the  equation  for  diffusion  with  drift  is 


du  du  <92u 

dt  C  dx  dx2 


(223) 


and  the  respective  Wigner  equation  of  motion  is 


dwu  dwu 

dt  +c  dx 


D  d2Wu 
2  dx2 


2  Dp2Wu 


(224) 


However,  no  generality  is  lost  by  taking  the  drift  term  equal  to  zero  because  if  u(x,  t )  solves  the  no 
drift  equation,  Eq.  (214),  then  u(x  —  ct,t)  will  solve  the  equation  with  drift.  Similarly,  if  Wu(x,p,t) 
satisfies  Eq.  (222)  then  Wu(x  —  ct,p,t)  satisfies  Eq.  (224). 


16.1  Green’s  function 


Schroeder  equation.  Suppose  we  want  to  solve  the  initial  value  problem  for  the  Schroeder  equation. 
That  is,  given  ip(x,  0)  we  want  ijj(x,t),  where  t  >  0.  The  solution  is 


ip(x,t)  =  /  G^(x,  x  ,  0)ip(x',  0)dx' 


(225) 


where  G^(x,  x\t)  is  the  Green’s  function, 


Gjj(x,  x',t)  =  exp 

V  4:iriat 


A  2 


(x  —  x') 
Aiat 


(226) 


In  momentum  space  the  initial  value  problem  becomes  particularly  easy.  From  Eq.  (217)  we  have 


(p{p,  t)=e  iap2t(j){p,  0) 


(227) 
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Now  consider  the  same  problem  for  the  Wigner  distribution,  that  is,  given  W(x,p,  0)  we  want 
W(x,p,t).  From  Eq.  (?????)  it  follows  that 

W^(x,p,t)  =  W^(x -2apt,p,0)  (228) 

a  result  first  obtained  by  Wigner  and  Moyal.  Thus  a  remarkable  simplification  is  achieved  in  phase 
space.  But  furthermore  in  phase  space  we  understand  what  is  going  on.  It  shows  that  as  time 
progresses  the  phase  space  point  moves  with  a  constant  velocity  in  the  x  direction  but  does  not  move 
at  all  in  the  p  direction.  The  velocity  in  the  x  direction  being  2 ap. 

Diffusion  equation.  Now  consider  the  diffusion  equation.  Using  the  Green’s  function  approach  one 

has 

u(x,t)  =  J  Gu(x,  x',  t)u(x' ,  0)dx>  (229) 

where 

G“(I’  x’’t)  =  vkm  “p  [~  <230) 

and  in  momentum  space 

U(p,t)  =  e~DP2tU(p,  0)  (231) 

Now  consider  the  Wigner  distribution.  One  can  show  that 

Wu(x,p,  t)  =  -j^=e^DpH  J  exp  ~^X  2pt  ^  Wu(x  ,p,())dx’  (232) 

In  Figs.  19-22  we  show  the  Wigner  distribution  computed  at  times  t  =  0.01,0.1,1,10  respectively, 
and  with  D  =  100. 
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Figure  19:  Wigner  distribution  of  the  Green’s  function  for  the  diffusion  equation,  for  t  =  0.01. 


Figure  20:  Wigner  distribution  of  the  Green’s  function  for  the  diffusion  equation,  for  t  =  0.1. 
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0.02  N 


Figure  21:  Wigner  distribution  of  the  Green’s  function  for  the  diffusion  equation,  for  t  =  1. 


Figure  22:  Wigner  distribution  of  the  Green’s  function  for  the  diffusion  equation,  for  t  =  10. 
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Thus  we  see  the  Wigner  distribution  shows  a  significant  physical  difference  between  the  two  Green’s 
functions  though  superficially  the  original  wave  equations  and  Green’s  functions  are  very  similar.  In 
the  case  of  Wq^(x,p,  t),  Eq.  (234)  shows  that  each  spatial  point  gets  transformed  by  a  translation 
phase  space,  the  translation  being  x'  — >  x  +  2 apt.  But  for  the  case  Gu(x,  x',  t )  we  see  that  each  point 
gets  spread/contracted.  The  position  spreads  but  the  momentum  contracts! 


Example.  We  take  a  specific  example  and  we  work  out  the  two  cases  side  by  side.  Take 


ip{x,  0)  = 


1 


(27TCT2)1/4 


exp 


u(x,  0)  =- 


1 


:  exp 


(x  -  x0f 
4  a2 

(x  -  xof 
2a2 


+  ipox 

(237) 

l2' 

(238) 

\Z2tut2 

and  we  note  that  the  densities  in  both  cases  are  the  same.  The  respective  initial  Wigner  distributions 
are  calculated  to  be 


Wip(x,p,  0)  =  -  exp 
7 r 


{x  -  x0y 


2a2 


-2a2(p-p0y 


Wu(x,p,  0)  = 


1 


exp 


0 X-XQ )2  2  2 


<7^ 


(J  p 


(239) 

(240) 


27Ty/7rcr 

Thus,  initially,  as  expected  both  Wigner  distributions  are  essentially  the  same.  However  as  time 
evolves 


1 


W^(x,p,t)  =  exp 

o— (2  Dt.+o2)p2 


Wu(x,p,t)  = 


(x  —  2  apt  —  xo)z 

2^ 

(x'  —  x)2 


A^aV2Di 

1  1 


exp 


-2a2(p-p0)2 

( x '  —  xq)2"1 


:  exp 


47t5/2  \j2Dt  +  a 2 

where  for  the  last  step  we  have  used 

f  e-(x(x-xi)2-P(x-x2)2  dx  =  r~  ^ 


2Dt  \ 
(x  -  x0)2 
2Dt  +  a2 


exp 


<7^ 


dx 


—  (2  Dt  +  (72)p 2 


;  exP 


a(3 

ex  j3 


(xi  -  x2y 


(241) 

(242) 

(243) 

(244) 


Oi  T  fd 

We  now  discuss  the  physical  meaning.  In  the  case  of  the  Schroeder  equations  the  Wigner  distributions 
is  just  rotating.  However,  for  the  diffusion  case  it  is  spreading  in  the  x  direction  and  contracting  in 
the  p  direction. 


17  Example:  diffusion:  §  +  eg  = 


We  show  the  derivation  of  the  equation  of  motion  for  the  Wigner  distribution  for  the  diffusion  equation. 
We  work  out  the  case  of  diffusion  with  drift, 


du  du  _  nd2u 
dt  C  dx  dx2 


(245) 
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where  u  =  u(x,  t )  is  the  field,  c  the  drift  coefficient,  and  D  the  diffusion  coefficient.  To  apply  our 
method  we  first  rewrite  the  equation  as 


d  d  „  d2  1  ,  , 

m  +  cai~Da^\  "(l’t)  =  0 

We  now  apply  our  method  and  obtain  two  equations  for  Z(x,p,t,uj) 


[At  +  cAx  -  DA2x]  Z(x,p,t,u)  =  0 
[Bt  +  cBx  -  DB%]  Z(x,p,t,u )  =  0 


(246) 


(247) 

(248) 


Expanding  the  operators  we  have 


1  <9  c  d  D  d2 

2d~t~lU+2di~lCP~^d^ 
Id  c  d  D  d2 

2M  +  iU,+  2di  +  iCp-  4ft? 


+  Dp 2  +  iDp 


d_ 

dx 


+  Dp 2 


iDp 


d_ 

dx 


Z(x,p,t,uj )  =  0 
Z(x,p,t,  uj)  =  0 


(249) 

(250) 


We  add  the  two  equations  in  order  to  have  a  real  equation  for  the  Wigner  Z(x,p,t,oj) 


d 


d  D  d2 


dt  C  dx  2  dx 2 


+  2  Dp2 


Z(x,p,  t,  lo)  =  0 


(251) 


This  is  the  equation  of  motion  for  Z(x,p,t,uj).  However,  since  u>  does  not  appear  in  the  equation,  we 
can  integrate  it  out  to  obtain  an  equation  for  the  standard  Wigner  distribution,  W(x,p,  t), 

dW  dW  D  d2W 


dt 


+  c 


dx  2  dx2 


-  2 DpzW 


(252) 


18  Example:  Heat  Equation 

We  consider  the  heat  equation 

d2u(x,t )  _  du(x,t) 
dx2  dt 

where  a  is  a  real  constant.  We  have  hence 

® AXZUU  —  AtZuu 
aBxZu,u  —  BtZuu 

Adding  and  subtracting  these  equations  we  have 

cl  [A2  ±  B2]  ZUjU  =  [At  ±  Bt]  Zu,u 

We  choose  the  plus  sign  to  obtain 


(253) 


(254) 

(255) 


(256) 


a  [Ax  +  B2  ZU:U  — 


dZ, 


u,u 


dt 


(257) 
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But 


(258) 


and  hence  we  have  that 


A*  +  Bl=12&+^ 


a  d  Zu  u  2  ry  9ZU  u 

r~  +  2aP  ZU,U  = 


2  dx 2  '  dt 

Since  to  does  not  appear  in  the  coefficients  we  can  integrate  it  out,  and  we  obtain 


ad2Wu,u  2  dWu 

+  2  ap  WU)U  = 


(259) 


(260) 


2  dx 2  '  8t 

An  interesting  issue  is  whether  such  an  equation  can  be  put  into  a  form  that  is  close  to  the  Liouville 
sin  equation  as  was  done  with  Schroedinger’s  equation.  This  is  currently  being  investigated. 


19  Example:  Burger’s  Equation 


For  the  Equation 


dx2 

The  same  derivation  as  above  leads  to 


d2u(x,t)  du(x,t)  du(x,t) 
a — — - 1 - ^ - 1-  c — - =  0 


dt 


dx 


■  d2Zu 


2  dx2 


—  +  2 ap  Zu  u  -  c- 


dZu_n,  dZ,„ 


dx 


dt 


20  Example:  Wave  Equation 

We  now  consider  the  ordinary  wave  equation 

'  d2  _  1_  8 ^ 
dx2  c 2  dt2 

Applying  the  methods  developed  we  have  that 


u(x,t)  =  f(x,t) 


1 


A2  -  ±A2 
c2 


1 


B  ^B't 
cz 


ZU:U(x,p,t,u)  =  Zfj(x,p,t,uj) 
Substituting  the  operators,  expanding  and  collecting  one  has 


LZUtU(x,p,t,uj)  =  Zfj(x,p,t,u) 


(261) 


(262) 


(263) 


(264) 


(265) 
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where 


L  =  — 
16 


d 4 

dx 4 


1  (  d4 


c2  1  dx2dt 2 


+ 


+ 


d4  \  1  94' 

dt.2dx2  )  c4  df4 

d 3  <93 


+  2^ 


4c2 

dx2  dt 

p  dxdt2  +pa 

a2 

1 

■  u;2 

a2 

p2  a2 

dx 2 

2 

dx 2 

~  Adt2 

<93 


d3 


—  LO- 


dt2dx  dtdx2 


a 2 


dxdt 


+  pco 


a2 


dtdx 

a2 


+  2 


P 


J  2  2 
2?  W 


21  Classical  Wave  Equation 


We  want  to  apply  the  same  method  developed  for  the  case  of  ordinary  differential  equations  to  the 
classic  wave  equation 

u(x,t)  =  f(x,t)  (266) 

One  can  prove  with  the  same  considerations  used  for  ordinary  differential  equations  that  the  equation 
for  the  four  dimensional  Wigner  distribution  is 


a2  _  i  a2 

dx 2  c 2  dt 2 


—A2 

cz 


r2  r2 
-Hr  2 
cz 


-K-ip,ipAi  pi  cu) 


Kf,f(x,pA,v) 


(267) 


However  we  believe  that  it  is  impossible  to  obtain  a  differential  equation  for  the  three  dimensional 
Wigner  distribution  in  this  case.  One  can  convince  oneself  of  this  by  attempting  to  do  so  directly  by 
the  same  methods  that  have  been  applied  to  the  Schroeder  equation.  Alternatively  one  attempt  to 
get  it  is  by  integrating  out  u j  from  Eq.  (267).  But  it  is  not  possible  to  integrate  out  co  to  obtain  an 
equation  for  W(x,p,t).  This  is  due  to  the  fact  that  the  operators  At  and  Bt  contain  oj. 


22  Random  systems 

Historically,  stochastic  processes  methods  have  been  developed  for  the  time-invariant  stationary  case. 
However,  typically  in  nature  stochastic  processes  are  nonstationary,  but  very  little  work  has  been 
done  on  this  case  because  of  the  difficulties  involved.  Nonstationary  processes  come  about  in  two 
general  ways.  First,  is  that  the  physical  parameters  of  the  process  can  vary  in  time  and  secondly  even 
for  a  process  that  is  generally  thought  to  be  stationary  there  is  a  nonstationary  part,  the  transient, 
which  is  usually  neglected,  but  which  is  very  important  and  interesting.  We  believe  that  it  is  usually 
neglected  because  proper  methods  have  not  been  developed  to  study  it.  In  what  follows  we  present 
a  new  approach  that  can  handle  nonstationary  stochastic  systems.  We  now  motivate  the  need  for 
this  development  by  examples.  Many  systems  are  described  by  differential  equations  that  have  a 
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driving  function  that  is  random  and  in  particular  the  driving  function  is  often  white  noise.  The 
output  of  such  systems  is  also  “noise”,  but  it  is  generally  colored  noise  and  may  be  time- varying. 
Such  differential  equations  are  stochastic  differential  equations  and  of  course  the  most  venerable  is  the 
original  equation  of  Langevin  to  describe  Brownian  motion,  and  also  the  standard  Wiener  process  is 
such  a  system.  When  the  output  of  such  differential  equations  is  stationary,  then  the  description  is 
given  by  the  power  spectrum.  However,  since  it  is  generally  the  case  that  the  output  is  nonstationary 
at  least  for  a  certain  interval  of  time,  a  generalization  of  the  power  spectrum  is  needed.  We  present  a 
method  that  allows  one  to  study  situations  where  indeed  the  output  is  nonstationary  and  also  colored. 
Using  the  methods  presented  we  will  study  some  important  cases  and  in  particular  we  will  derive  new 
time- varying  spectral  properties  of  nonstationary  random  processes. 

Suppose  we  have  a  differential  equation  that  governs  the  process  x(t)A 


dnx(t )  dn  1x(t )  dx(t) 

dn i -  1 i - i —  *  *  *  d\ +  dQX{t ) 

dtn  dtn~x  dt 


m 


(268) 


and  where  fit)  is  white  Gaussian  noise  with 


E[f(t)}  =  0  (269) 

Rf(h,t2)  =  £[/(*i)/(t2)]  =  N0S(ti  -  t2)  (270) 

and  where  E[  ]  is  the  ensemble  averaging  operator.  We  consider  systems  defined  by  Eq.  (268)  where 
we  take  the  coefficients  to  be  constants.  Our  aim  is  to  understand  the  spectrum  of  the  output,  x(t). 
As  mentioned,  the  spectrum  will  in  general  be  colored  and  nonstationary.  We  use  the  Wigner-Ville 
approach  for  describing  such  time-varying  spectra.  The  Wigner  distribution  of  a  deterministic  signal 
x(t)  is  [?,  ?] 

Wx{t,Lu)  =  —  [  x*(t  —  r/2)x(t  +  T/2)e~lTbJdT  (271) 

27 T  J 

The  Wigner-Ville  spectrum  of  a  random  process  x(t)  is  defined  as  the  ensemble  average  of  the  Wigner 
distribution  of  x(t)  5  [?,  ?] 

Wx(t,uj)  =  ^  j  E[x*  {t  -  t /2)x{t  +  t /2)\e~lTU dr  (272) 


23  Nonstationary  stochastic  system 


Our  aim  is  to  obtain  the  equation  of  motion  for  the  Wigner  spectrum  for  a  stochastic  process  that  is 
the  solution  of  an  ordinary  stochastic  differential  equation 


dnx(t )  dn  1x{t)  dx(t) 

dry : -  +  dfi —  i : - ; —  •  •  •  H-  d\ +  dnx(t) 

dtn  dt71-1  dt  w 


fit) 


4We  use  bold  symbols  for  stochastic  quantities. 
5We  use  the  following  simplified  notation  f  = 


(273) 
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where  f(t)  is  a  given  stochastic  process.  One  could  solve  this  equation  for  x(t)  and  then  find  the 
ensemble  average  as  needed  to  calculate  the  Wigner  spectrum  in  Eq.  (272),  and  then  do  the  integration 
as  required.  This  is  generally  a  difficult  procedure.  We  have  developed  a  method  that  is  more  direct: 
we  obtain  the  equation  of  evolution  for  the  Wigner  spectrum  and  solve  it  [?,  ?,  ?,  ?].  To  accomplish 
that  we  write  the  differential  equation  in  polynomial  form 


P(D)x(t)  =  f(t) 


(274) 


where,  as  usual,  D  and  P(D )  are  respectively 

«  =  4 

dt 

P(D )  =  anDn  +  an_iT'l_1...  +  m  D  +  a0 

The  differential  equation  for  the  Wigner  spectrum  Wx(t,Lo),  of  x(t),  is  then  given  by 

P*(A)P(B)Wx(t,u)  = 


where 


Id.  Id. 

A=29tr,U’  '  B=2»+JU’ 


For  the  case  of  time  dependent  coefficients 


,  dnx(t )  .dn  1x(t)  dx(t ) 

+an-i(t)  x  +  ••  .  +  ai(f)— —  +a0{t)x(t)  =  f(t ) 


(275) 

(276) 

(277) 

(278) 

(279) 


dtn  —v-'  dtn- 1  dt 

where  ao(t), . . .  ,an(t )  are  time-varying  deterministic  coefficients.  We  rewrite  Eq.  (279)  in  polynomial 
notation 

P(D,t)x(t)  =  x(t)  (280) 

where  now 

P{D ,  t)  =  an(t)Dn  +  an-i{t)Dn  +  . . .  +  a\{t)D  +  ao(t)  (281) 

The  equation  for  the  Wigner  spectrum  is 


P*{A,  S)P(B,  P)Wx(t,  u)  =  W  f(t,u) 
Eq.  (282)  is  a  partial  differential  equation  in  time  and  frequency. 


(282) 


24  Example:  The  Nonstationary  Wiener  Process 

We  now  apply  our  method  to  the  Wiener  process  [?] 

^  -  /<*> 


(283) 
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with  initial  conditions  *(0)  =  xq,  where  xq  is  a  stochastic  variable.  We  will  turn  the  system  on  at  a 
finite  time  and  hence  we  will  obtain  the  complete  solution  of  the  problem.  That  is  we  will  get  the  total 
solution,  including  the  transient  behavior.  First  we  point  out  that  for  Eq.  (283)  the  power  spectrum 
is  given  by 

Px{u)  =  ^L\  (284) 

V27T 

It  is  important  to  appreciate  that  the  1/ui2  spectrum  for  the  Wiener  process  is  reached  in  either  of 
two  different  ways: 

1.  The  initial  condition  is  set  to  zero  at  t  =  0  and  by  taking  t  =  oo  in  the  general  solution  of  Eq. 
(283); 

2.  Alternatively  by  setting  the  initial  condition  to  zero  at  t  =  — oo  and  this  has  the  effect  of  the 
passage  of  an  infinite  amount  of  time. 


For  both  types  of  initial  conditions  the  stationary  spectrum  is  reached  because  an  infinite  amount 
of  time  has  passed.  Using  our  method  we  will  obtain  the  full  solution,  that  is  the  full  evolution  of  the 
instantaneous  spectrum.  Rewriting  Eq.  (283)  in  the  polynomial  form  of  Eq.  (274) 

Dx(t)  =  /(f)  (285) 


and  using  Eq.  (277)  we  obtain 

ABWx(t,  uj)  =  Wf(t,u ) 

Using  the  definitions  of  the  operators  we  have  that 


and  hence 


ld2WAt’U)  +ulWAt,u)  =  Wf(t,u) 


(286) 


(287) 


(288) 


4  dt’2 

where  Wf(t,u>)  is  the  Wigner  spectrum  of  f(t).  If  /(f)  is  white  Gaussian  noise  with  autocorrelation 


then  one  can  readily  show  that 


and  hence 


Rf(T)  =  NqS(t) 


wf (*■">  =  T, 


ld2W„(t,u)  2TI7  u  ,  JV„ 

I  W  +W  W^=2i 


The  exact  solution  to  Eq.  (291)  is 


Wx(t,  to)  =  [i  _  cos  2 ut\  +  E[xo\  gjn  2 cot,  t  >  0 

27ruU  7 TOJ 


(289) 

(290) 

(291) 

(292) 
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Figure  23:  The  instantaneous  spectrum  of  the  Wiener  process  defined  in  Eq.  (283).  The  spectrum 
grows  to  infinity  at  lu  =  0  and  reaches  the  singularity  in  the  classical  frequency  spectrum  that  has  a 
1/a ;2  distribution.  Notice  that  the  Wigner  spectrum,  given  in  Eq.  (293),  is  nonnegative. 

and  Wx(t,u)  =  0  when  t  <  0.  For  E[x q]  =  0  one  has 

Trr  \  N°  n  o  u  No  sin2  tot 

Wx(t,u )  =  - - ^  1  -  cos  2a ;t\  = - 5— 

2hujz  7T  ojz 

In  Fig.  23  we  show  the  transient  of  the  Wiener  process. 

It  is  interesting  to  note  that  for  10  — >  0 

lim  Wx(t,uj)  =  lim  — ^  [1  —  cos2a>f]  =  — -t2 

a;— >0  ’  u.^0  2-KUJ2  L  J  7 T 

and  hence  at  lu  =  0  the  instantaneous  spectrum  grows  to  infinity  and  approaches 
the  classical  power  spectrum  of  Eq.  (284)  with  a  second  order  in  time. 

24.1  Direct  Solution 

We  now  show  how  to  obtain  the  same  solution  by  the  direct  method  which  is  essentially  a  brute  force 
method.  The  Wiener  equation  is  first  rewritten  as 

=  u(t)f(t )  (295) 


(293) 


(294) 

the  singularity  in 
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where  u(t)  is  the  Heaviside  step  function,  defined  as 


u(t)  = 


0,  t  <  0 
1  t  >  0 


The  use  of  the  step  function  makes  the  calculation  much  easier.  Now,  we  write  the  solution  to  Eq. 


(295)  as 


x(t)  =  u(t)  xq  +  f  u(t')f{t')dt' 

L  Jo 


To  obtain  the  Wigner  spectrum  we  will  make  use  of  the  following  relation  of  the  Wigner  spectrum 
with  the  autocorrelation  function,  Rx(ti,t2), 

Wx(t,  to)  =  J  Rx(t  +  r/2,  t  -  t /2)e~iTUJdr  (298) 

Hence,  we  need  to  evaluate  the  autocorrelation  directly 


Rx(h,h)  =  E[x(h)x*  (t2)]  (299) 

pt  i  rt2 

=  u{t\)u{t2)  E[x o]  +  /  /  Rf(t\ ,  t2)dt\dt'2  (300) 

L  Jo  Jo 

where  the  cross-terms  have  disappeared  because  E[f(t)]  =  0.  By  substituting  R.fit^,  t'2)  from  Eq. 
(270)  one  obtains 

Rx(h,t 2)  =  u(ti)u(t2)  mm(ti,t2)  (301) 

Now,  one  applies  Eq.  (298)  to  obtain 

Wx(t,uj)  =  j  u{t  +  r/2)u(t  —  r/2)  (min(t  +  r/2,  t  —  r/2)  +  Efccg])  e~lT^dr  (302) 

and  after  some  manipulation  and  simplification  (shown  in  the  next  Section)  one  has 

Wx(t,  to)  =  9  [1  —  cos  2c ut]  +  sjn  2 uit,  t  >  0  (303) 

2tTUJZ  7 TiO 

which  is  precisely  Eq.  (292).  We  note  that  Wx(t,uj)  =  0  for  negative  times. 

24.2  Direct  derivation  of  the  Wiener  process 

We  now  give  the  details  of  the  derivation  of  the  direct  approach  leading  to  Eq.  (272).  We  start  from 
Eq.  (299).  By  substituting  Eq.  (270)  into  Eq.  (300)  we  have 

Rx(h,t2)  =  u(ti)u(t2)  E[ Xq\+  [  [  N08(t[  -  t^dt^dt^  (304) 

Jo  Jo 

Now,  if  t\  <  t2  we  have  that 


Rx(ti,t2)  =  u(ti)u(t2)  [E[x 0]  +  N0ti] 
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while  if  f  2  <  f  1 

Rx(h,t2)  =  u(ti)u(t2)  [E[x q]  +  N0t2] 
which  can  be  combined  into  one  expression, 

Rx{ti,t2)  =  u(ti)u(t2)  [E[xq]  +  A^min(ti,t2)] 

Also,  it  is  convenient  to  write, 

Rx(h,t2 )  =  u(ti)u(t2)  [E[x q]  +  No  (u(ti  -  t2)t2  +  u(t2  -  fi)ti)] 
Now  we  apply  Eq.  (298)  that  links  the  autocorrelation  to  the  Wigner  spectrum 


Wx(t,u)  = 


E[x20 ] 

2t r 


u(t  +  t /2)u{t  -  t /2)e~iru} dr  + 


h 


— -  [  u(t  +  t /2)u{t  —  t /2)u(r){t  —  t /2)e  lTLJdr  + 
27 t  J 

" - v - ' 

h 

— -  [  u(t  +  r/2)u(t  —  T/2)u(r)(t  +  t /2)e~lTU} dr 
27 t  J 

' - - - - - ' 

h 


We  solve  the  three  integrals  separately.  Starting  with  the  first,  we  have  that 

h  =  u{t)- 


E[x2]  r2t 


— /  e~iruJdT 

2^ r  J-2t 


,  .  „r  9l  sin  2 t,io 

=  u(t)E[x  5]- 


7 TLO 


For  the  second 


N0 


f2 1 


h  =u(t)-E  (t-T/2)e~'™dT 

2tt  Jo 


which  after  a  few  calculations  becomes 


h  =  »(*)N 


t- 


VjJ 


1 

—i2tu  1 

(  i 

2  A 

1  ' 

2 

e  1 

iuj  J 

_  A 

With  a  similar  procedure  one  has  that 


h  =  u(t )^°  /°  {t  +  T/2)e-™dT 
=  u{t)^-j\t-T/2)e-dr 


which  evaluates  to 


/3  =  u(i) 


Ao 

2tt 


1  -  e~i2tuJ  1 


707 


V-)-  4 


(306) 

(307) 

(308) 

(309) 

(310) 

(311) 

(312) 

(313) 

(314) 

(315) 

(316) 

(317) 

(318) 
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Note  that 


jVn 

h  +  h  =  u(t )  [1  -  cos  2fw] 


(319) 


and  now  adding  to  I\  we  obtain 

Wx(t,ijj)  =  u(t) 


N0  ,  , 

- 7T  1  —  cos  2cuf  H - —  sin  2 uit 

2itujz  TTU) 


(320) 


which  is  Eq.  (303)  of  the  text.  We  note  that  the  phase  space  method  used  before  is  much  easier. 


25  Example:  The  full  exact  solution  to  the  Langevin  Equation 

Our  aim  is  to  show  how  to  find  the  transient  spectrum  of  the  Langevin  equation 

dv(t) 


dt 


+  pv(t)  =  /(f) 


(321) 


where  /(f )  is  white  Gaussian  noise.  We  apply  our  method  by  first  rewriting  the  Langevin  equation  in 
polynomial  form 

[D  +  0]v(t)  =  f(t)  (322) 


Using  our  method  the  equation  for  the  Wigner  spectrum  is  now 

[A  +  P}[B  +  P]Wx(t ,  u)  =  Wf(t,  lo) 

We  substitute  for  the  operators  A  and  B  from  Eq.  (278)  to  obtain 

\  w)  +  P§jW9{t,  u)  +  (P2  +  U2)Wx(t,  U>)  =  Wf(t,  (j) 


(323) 


(324) 


where  Wf(t,u> )  is  the  Wigner  spectrum  of  the  input  random  process  /(f). 

The  solution  to  Eq.  (324)  can  be  found  with  the  standard  methods  for  differential  equations.  It  is 


N0  1 


N0  e 


—2/3* 


(cos 2ut  —  lu/P sin 2o/f)  ;f>  0  (325) 


27T  P2  +  LO2  27 r  P2  +  LU2 

and  where  Wv(t,u)  =  0  for  negative  times.  In  the  solution  we  are  considering  a  random  initial 
condition  Vq,  that  is 

v(0)  =  v0  (326) 

This  solution  allows  us  to  understand  the  transient  behavior  of  the  Langevin  equation.  In  Fig. 
24  we  show  a  typical  case.  We  see  that  the  instantaneous  spectrum  has  an  overshoot  effect  at  the 
beginning  of  the  transient.  Also  the  spread  of  the  solution  is  larger  for  small  times.  Then  the 
solution  rapidly  approaches  the  stationary  spectrum,  the  well  known  Lorentzian  distribution.  The 


54 


t(s) 


Figure  24:  The  complete  instantaneous  spectrum  of  Brownian  motion,  given  in  Eq.  (325).  Notice  the 
overshoot  at  the  beginning  of  the  transient,  and  the  increased  spread  of  the  Wigner  spectrum.  As 
t  — >  oo  the  Wigner  spectrum  reaches  the  stationary  solution,  that  is  the  Lorentzian  distribution  of 
Eq.  (327). 
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behavior  of  the  frequencies  in  the  transient  suggests  the  possibility  of  designing  systems  so  to  have 
smooth  transitions  of  the  instantaneous  spectrum  from  zero  to  the  stationary  solution.  This  can  be 
a  fundamental  issue  in  minimizing  the  electromagnetic  interference  of  an  electronic  device  that  is 
suddenly  turned  off,  or  in  reducing  the  risk  of  breaking  down  an  engine  or  a  vibrating  structure  by  a 
sudden  rough  start.  These  two  physical  situations  are  just  an  example  of  the  importance  of  transients 
in  designing  and  understanding  a  physical  systems.  We  note  that  as  time  goes  to  infinity  we  have  that 


lim  Wv(t,  oj) 
t->  o 

which  is  the  classical  Wang-Uhlenbeck  result. 


Nq 

Po  +  <^2 


(327) 


26  Example:  Quantum  Langevin  Equation 


The  quantum  Langevin  equation  is 

m^+m^  +  r(,)  =  m  (328) 

where  x  is  the  position  operator,  U(x),  the  external  potential,  and  £(f)  is  the  noise  operator  that 
satisfies 

([£(*)>  £(*')]+)  =  ^  J  we!a,(^'')  coth  ^div  (329) 

If  one  considers  Eq.  (328)  as  a  classical  type  equation  then  it  is  called  the  quasi-classical  Langevin 
equation.  In  such  a  case  one  replaces  the  operators  by  ordinary  variables 

rn<i  ^  +  v ,(*)  =  £(*)  (33°) 
and  the  autocorrelation  function  is  then 

R(t)  =  —  [  u>elUT  coth  -^-dio  (331) 

V  '  ^  i-oo  2  kT  V  ; 

We  consider  here  the  case  of  no  external  potential  and  rewrite  Eq.  (328)  as 


fi(t)+Pp(t)  =  €(t) 


(332) 


with 


and 


R^t)  =  2 DZ  j 


ujelu,T  coth  Zujdu) 


Z  = 


h 

2 kT 


Using  the  standard  Wiener-Khinchen  theorem  the  power  spectrum  is  given  by 


(333) 

(334) 


S^(io)  =  2DZlo  coth  Zlu 


(335) 
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We  note  that  as  Z  — >  0 


(336) 


|ni  i??(r)  =  2 D5{t) 

Using  our  general  approach  we  have  that  the  differential  equation  for  the  Wigner  spectrum  is 


n— 

4  dt 2  dt. 


DZ 


7^  +  +  /?2  +  w2  |  Wp(t,  u)  —  — — wcoth  Zlo 


The  exact  solution  is 


Wp(t,w)  = 


f32  +  cj2  L 

1  DZ 


_  2 7, 

1  —  e  cos  2u;t 


w  coth  Za; 


_27f 

1  —  e  "*  cos  2wt 


(337) 

(338) 

(339) 


P2  +  U2  77 

We  see  that  Eq.  (339)  can  now  be  thought  of  as  the  generalization  of  the  Wang  and  Uhlenbeck  process 
for  quantum  noise. 

The  quantum  Wiener  process  is  obtained  by  taking  7  =  0  and  hence  the  differential  equation  is 


p{t)  =  £(*) 


(340) 


The  governing  equation  is 


1  d2Wv(t,L0)  2tt7  1  x 

4 — ^r^+^2M rv(tM 


DZ 


-lo  coth  Zlo 


77 


(341) 


giving 


Wp(t,u) 


DZ 

- oj  coth  Zu>(  1  —  cos  2 ujt) 

77 


2  DZ 


(jo  coth  Zu  sin2  out 


77 


(342) 

(343) 


27  Example:  Time-variant  random  systems 

In  many  physical  situations  a  random  process  can  be  seen  as  the  output  of  a  random  system  whose 
parameters  change  with  time.  The  variation  can  be  due  to  a  sudden  breakdown,  aging,  or  cyclic  vari¬ 
ations  of  temperature,  humidity,  and  other  physical  quantities  that  can  influence  the  system  behavior. 
In  these  cases  the  process  is  inherently  nonstationary. 

We  now  consider  the  case  addressed  by  a  time-dependent  Brownian  motion  process 

+  0  (t)v(t)  =  f(t )  (344) 

We  rewrite  it  in  polynomial  notation 

[D  +  P(t)]v(t)  =  f(t )  (345) 
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where  we  take  (3(t)  to  be  time-dependent  with  a  linear  law 


Pit)  =  Po  +  et 


(346) 


and  also  we  take  e<  1.  These  conditions  imply  a  slow  variation  with  time.  We  now  transform  the 
equation  in  time  to  an  equation  in  time-frequency,  using  Eq.  (282)  to  obtain 


1  d 2 


+ 


d 2 


d 


d 


+  P(t)  wt  +  eu;  W - ^  P  (t)  +  u 


Wv(t,u )  = 


4  <9f2  4  dw2  '  r'x"/dt  '  ~  dtu 

Taking  into  account  that  e  is  small,  we  neglect  the  terms  containing  e, 

7A17  +  Pi^^l  +/?2(t)  +^2  kEv(t,w)  =  ^ 


iVo 

2tt 


(347) 


(348) 


_4df2  ' 

Since  there  are  no  derivatives  with  respect  to  a;,  we  can  solve  this  equation  as  an  ordinary  differential 
equation  with  respect  to  time.  We  have  to  choose  the  initial  conditions,  which  we  fix  using  the 
following  reasoning.  We  assume  that  the  system  was  started  at  — oo  with  (3  =  Po  and  hence  it  has 
reached  equilibrium.  At  t  =  0  we  turn  on  the  time  dependence  as  given  by  Eq.  (346).  Therefore  the 
initial  conditions  to  Eq.  (348)  are 

Ao 


Wv{  0,cu)  = 


Pn+U2 


^Wv(t,u)  =  0,  t  =  0 


(349) 

(350) 


The  first  initial  condition,  Eq.  (349),  is  the  Wigner  spectrum  of  the  stationary  solution.  This 
stationary  distribution  has  been  reached  because  an  infinite  amount  of  time  has  passed  since  the 
system  was  turned  on  at  t  =  — oo.  As  a  consequence  of  being  in  a  stationary  phase  before  t  =  0  there 
must  be  no  change  in  time  of  the  instantaneous  spectrum.  This  condition  is  stated  by  Eq.  (350). 

Using  these  initial  conditions  Eq.  (348)  can  be  easily  solved  with  numerical  algorithms.  In  Fig. 
25  we  plot  the  solution  Wv{t,  uj).  In  Fig.  26  a  numerical  simulation  of  the  Wigner  spectrum  is  shown. 
The  agreement  confirms  our  approach. 

The  behavior  of  the  instantaneous  spectrum  is  in  accordance  with  intuition.  The  shape  of  the 
Wigner  spectrum  is  in  fact  very  similar  to  a  Lorentzian  distribution,  for  any  given  time.  A  similarity 
with  the  lowpass  filter  represented  by  the  RC  circuit  can  be  used  to  explain  the  change  in  spread  that 
the  instantaneous  spectrum  exhibits.  An  RC  circuit  driven  by  white  noise  has  the  same  equation  of 
Brownian  motion.  In  this  case  we  can  think  of  a  time- varying  RC  circuit,  where  the  time  constant 
r  =  RC  is  a  function  of  time,  r  =  r(t).  It  is  well  known  that  when  r  is  large,  the  RC  filter  has  a  strong 
lowpass  behavior,  which  corresponds  to  a  small  spread  of  the  Lorentzian  spectrum.  On  the  contrary, 
when  r  is  small,  the  RC  circuit  has  a  weak  lowpass  behavior.  Since  the  time  constant  is  related  to  the 
/3  coefficient  by  an  inverse  proportionality,  the  increase  of  the  /?  coefficient,  which  we  assumed  for  our 
case,  turns  out  to  be  a  decrease  of  the  r  constant.  Therefore  we  expect  the  time- varying  system  under 
study  to  become  a  weaker  lowpass  filter  as  time  goes  by,  which  is  precisely  what  happens  in  Fig.  25. 
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Figure  26:  Simulated  Wigner  spectrum  of  a  time-varying  Brownian  motion. 


28  Example:  oscillator  with  constant  coefficients 


We  consider  the  case  of  a  harmonic  oscillator  with  constant  coefficients  and  random  Gaussian  noise 
as  input, 

d2x  dx  9  „  , 

^+2„- +  „„*  =  /(()  (351) 

where  we  choose  a  Gaussian  noise  /(f)  with  zero  mean  and  autocorrelation  Rf(r )  =  Nq5(t).  Our 
intention  is  to  study  how  the  frequency  in  the  random  output  x(t)  of  the  systems  are  distributed  in 
time,  that  is  the  instantaneous  power  spectrum  as  defined  by  the  Wigner  spectrum.  Of  course  for 
this  case  we  expect  the  spectrum  to  be  constant  in  time.  First  we  consider  the  problem  by  standard 
methods  and  then  apply  our  method  and  compare.  Subsequently  we  consider  the  problem  where  we 
make  one  of  the  constants  time  dependent,  and  of  course  the  standard  method  can  not  deal  with  that 
situation. 

28.1  Standard  result 

We  recall  the  standard  result  for  the  classical  power  spectrum  Gx{u)  of  x(t),  defined  as  usual  as  the 
Fourier  transform  of  the  autocorrelation  function  Rx(t) 

Gx{u)  =  \H(u)\2Gf{u)  (352) 

where  H{uj)  is  the  transfer  function  of  the  system,  simply  obtained  by  evaluating  the  polynomial 
version  of  Eq.  (268)  in  iu,  that  is 


H(  u)  =  P(iu)  (353) 

(Notice  that  with  constant  coefficients  the  polynomial  is  not  a  function  of  time  anymore).  Since  the 
power  spectrum  of  the  Gaussian  noise  is  constant  and  equal  to  Gf(uj)  =  No,  we  obtain  from  Eqs. 
(352)-(353) 


Gx{u)  = 


N0 


(uJq  —  to2)2  +  4/i2cu2 


(354) 


This  is  the  result  obtained  by  Wang  and  Uhlenbeck  [?].  In  Fig.  27  we  represent  the  obtained  power 
spectrum  Gx( u)  when  the  following  set  of  parameters  is  chosen 


l-l  =  1 

(355) 

UJO  =  27T 

(356) 

No  =  1 

(357) 
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Figure  27:  Power  spectrum  of  the  state  variable  x(t)  of  the  harmonic  oscillator  of  Eq.  351.  Notice 
the  bandpass  behavior  of  the  function. 


We  point  out  the  bandpass  behavior  of  the  power  spectrum,  that  selectively  filters  the  input  frequencies 
contained  in  the  random  noise  f(t). 

28.2  Wigner  spectrum 

We  now  obtain  the  equation  for  the  Wigner  spectrum  by  applying  our  method.  We  rewrite  Eq.  (351) 
in  polynomial  notation 


[ D 2  +  2 /jlD  +  cuq]  x(t)  =  f(t) 

and  we  apply  the  transformation  to  the  Wigner  spectrum  domain,  having 


(358) 


[A2  +  2j.iA  +  Wq]  \_B~  +  2/i.B  +  Wq]  W x(t,  to)  —  W  f(t ,  u)  (359) 

This  is  a  partial  differential  equation  of  fourth  order,  that  can  be  written  in  the  following  general  way 


\  d4  d3  ,  d2  ,  d  ,  ' 

hdB  +  b3W  +  b2w+blm  +  bo 


Wx(t,u)  =  Wf(t,u ) 


(360) 


One  can  show  that  the  coefficient  bo  is  inversely  proportional  to  the  square  modulus  of  the  transfer 
function 


b0=(\H(u)\- 


-i 


(361) 
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This  means  that  in  general  Eq.  (360)  can  be  written  as 


<94  d3  d2  d 

hm 4 +  +  b2w+blm  +  ( ,F(w)r 


-1 


wx(t,uj)  =  Wf(t,u)  (362) 

As  can  be  easily  seen,  the  Wigner  spectrum  of  white  noise  is  a  (two-dimensional)  constant  function 


W  =  No  (363) 

We  then  use  this  result  in  Eq.  (362)  and  we  solve  the  equation  as  prescribed  by  our  method,  that 
means  we  look  for  the  solution  that  is  the  convolution  of  the  Green’s  function  with  the  input  function. 
This  solution  is  seen  to  be  constant  with  respect  to  time,  and  in  particular 

W9(t,u)  =  N0\H(u)\2  (364) 

Surprisingly,  we  obtain  the  same  solution  as  in  the  case  of  the  classical  power  spectrum.  But  there 
is  the  important  difference  that  here  we  are  evaluating  a  function  of  time  and  frequency.  The  result 
obtained  means  that  the  instantaneous  power  spectrum  of  the  solution  x(t)  is  constant  in  time.  A 
simple  and  intuitive  explanation  of  this  conclusion  can  be  given  by  considering  Eq.  (351)  as  a  filtering 
problem.  In  this  perspective  what  the  system  does  is  filter  the  input  noise  f(t)  with  a  bandpass 
filter  to  produce  an  output  signal  x(t).  But  since  both  the  input  random  noise  and  the  filter/system 
are  stationary  in  time,  then  also  the  output  is  stationary.  This  stationary  is  contained  in  the  time- 
invariant  property  of  the  instantaneous  power  spectrum  Wx(t,uS).  Also  it  can  be  easily  noticed  that 
the  property  obtained  in  Eq.  (364)  is  valid  for  any  system  defined  by  Eq.  (268). 


29  Example:  Harmonic  oscillator  with  time  dependent  coefficients 

We  now  consider  a  time-variant  harmonic  oscillator  described  by  the  following  differential  equation 

dznr  Hnr 

—  +  2»—  +  K(t)x  =  f(t)  (365) 

where 

K(t)  =  cd'o  +  et,  e<l  (366) 

This  is  a  slightly  perturbed  version  of  the  standard  harmonic  oscillator  that  shows  though  a  basic 
time-varying  behavior.  Our  intention  is  to  study  the  system  in  the  time  interval  0  <  t  <  T,  where  T 
is  large  enough  to  show  a  variation  in  the  K(t)  coefficient.  In  particular  we  choose 

K(T)  =  4A'(0)  =  4  ul  (367) 
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The  equation  for  the  Wigner  spectrum  is 


[A2  +  2 nA  +  K(E)\  [ B 2  +  2fiB  +  K(F )]  Wx(t,u)  =  Wf(t,u )  (368) 

This  equation  is  a  partial  differential  equation  with  varying  coefficients,  that  cannot  be  solved  as  an 
ordinary  differential  equation  (there  are  derivatives  with  respect  to  lo).  But  if  we  take  into  account 
the  fact  that  K(t)  is  slowly  varying,  we  can  approximate  it  by 


16  dt 4  2  dt3 


+  +  2  (*)  + 


dt 2 


+2 n  ( K(t )  +  u)2)  +  (/i  (t)  -  J2)2  +  4:fi2L02 


Wa 


Wf  (369) 


This  is  again  a  partial  differential  equation  but  contains  no  derivatives  with  respect  to  oj  and  can  hence 
be  solved  as  an  ordinary  differential  equation.  The  only  other  issue  is  to  set  proper  initial  conditions 
at  t  =  0.  Since  the  system  is  slowly  time-varying,  due  to  the  fact  that  e  is  very  small,  we  suppose 
that  its  instantaneous  spectrum  will  be  slowly  varying.  If  this  is  the  case,  then  the  instantaneous 
spectrum  at  t  =  0,  when  K(t  =  0)  =  lUq,  will  not  be  very  different  from  the  instantaneous  spectrum 
of  the  stationary  system  that  has  always  K(t)  =  cuq  (the  one  analyzed  in  Sect.  28).  Following  this 
approximation  we  choose 


Wx(  0,uj) 


dWx 

dt 


(0,w) 


N0 


Wn  - 


4  fi2LU2 


d2w*<*  ^_&wX/n  , 

dt 2  (0,w)  dt 3  (0,w) 


=  o 


(370) 

(371) 


The  solution  of  Eq.  (369)  with  the  set  of  initial  conditions  given  in  Eqs.  (370)-(371)  can  be 
obtained  analytically  using  a  power  series  expansion.  In  Fig.  28  we  represent  this  solution,  that  has 
been  instead  obtained  with  a  standard  numerical  integration  method.  This  is  necessary  since  the 
domain  of  integration  is  large  and  the  power  series  solution  should  be  computed  for  too  many  terms, 
far  beyond  the  machine  precision.  (The  problem  can  be  handled  with  a  symbolic  program  but  here 
we  use  a  numerical  integration  scheme  because  of  its  easy  and  robust  implementation). 

It  is  very  important  to  notice  that  the  computed  Wigner  spectrum  Wx(t,uj)  looks  exactly  as  we 
expected.  It  is  basically  a  bandpass  spectrum  as  in  Fig.  352  with  a  time-varying  central  frequency. 
Also  the  amplitude  of  the  spectrum  changes  in  time,  and  the  reason  is  that  this  is  embedded  in  the 
equation  defining  the  time- varying  harmonic  oscillator,  Eq.  (365).  The  smooth  transition  between  the 
initial  spectrum  at  t  =  0  and  the  final  spectrum  at  t  =  T  is  a  consequence  of  the  small  value  of  the  e 
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Figure  28:  Wigner  spectrum  of  the  solution  x(t)  to  Eq.  365. 

parameter,  that  makes  K (t)  =  Uq  +  et  a  slowly  varying  function.  Now  the  crucial  step  is  to  compare 
this  result  that  satisfies  our  intuition  with  the  estimated  Wigner  spectrum  obtained  by  simulation. 
This  is  done  in  Sect.  29.1. 

29.1  Comparison  with  simulations 

To  check  the  validity  of  the  solution  shown  in  Fig.  28,  we  compare  it  with  the  estimated  time-varying 
spectrum  obtained  by  simulation.  We  first  solve  Eq.  (365)  with  a  scheme  that  is  valid  for  a  stochastic 
differential  equation  (e.g.  an  Euler  scheme),  and  we  thus  obtain  an  approximation  to  x(t).  Then  we 
estimate  the  instantaneous  power  spectrum  Px(t,oj)  with  a  sliding  estimator  method.  This  method 
consists  in  truncating  the  signal  with  a  moving  window,  and  then  estimating  the  power  spectrum  at 
that  time  with  a  standard  technique.  Here  a  parametric  ARMA  estimator  has  been  used,  since  we 
know  that  the  truncated  signal  is  coming  from  a  second  order  differential  equation.  In  Fig.  29  we 
show  the  result  of  such  estimation.  It  can  be  noticed  that  our  approximation,  shown  in  Fig.  28  and 
the  estimated  spectrum  are  very  close  to  each  other.  In  particular  Fig.  29  shows  again  the  bandpass 
behavior  that  we  were  expecting  from  an  intuitive  point  of  view  and  that  our  approximation  confirmed. 
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Figure  29:  Numerical  estimation  of  the  Wigner  spectrum  of  the  solution  to  Eq.  (365).  Notice  the  qual¬ 
ity  of  the  approximation  obtained  in  Fig.  28  that  fits  in  an  excellent  way  the  estimated  instantaneous 
spectrum. 
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30  Clouds 


The  study  of  cloud  coverings  is  important  for  many  reasons,  among  them  is  that  we  want  to  remove 
the  clouds  from  images.  Our  aim  is  to  understand  certain  aspects  of  clouds  in  regard  to  their  statistical 
properties.  Considering  a  cloud  as  a  signal  we  ask  where  in  the  Fourier  transform  is  the  cloud-like 
information  stored.  In  particular,  is  it  in  the  spectral  amplitude  or  spectral  phase?  The  reason  this  is 
important  is  that  by  understanding  this  issue  it  may  help  in  designing  methods  for  the  denoising  of 
clouds,  that  is,  methods  that  eliminate  cloud  coverings  [?,  ?].  The  idea  of  studying  the  relative  impor¬ 
tance  of  spectral  amplitude  and  phase  of  a  signal  originated  with  the  classical  paper  by  Oppenheim 
and  Lim  [?]  and  subsequently  many  studies  have  been  made. 

We  will  start  by  studying  these  issues  on  artificially  generated  clouds. 


30.1  Generation  of  clouds 

For  a  two  dimensional  image,  s(x,y),  we  write  its  Fourier  spectrum  as 

F(ux,  u)y)  =  j  J  s(x,  y)e~lxuJx~iyuJydxdy  (372) 

and  consequently  the  power  spectrum  is  given  by 

PS((jJX,LOy)  =  \F(ux,ujy)\2  (373) 


One  way  of  generating  cloud-like  images  is  by  using  the  fact  that  their  Fourier  power  spectrum  is 
of  the  type 


Ps  i  i  M  y  )  — 


1 


,27 


(374) 


where  7  is  a  parameter  that  controls  the  cloud-like  appearance  of  the  image.  This  type  of  power 
spectra  is  observed  in  natural  clouds  [?,  ?,  ?,  ?].  When  one  wants  to  generate  artificial  clouds  it  is 
hence  reasonable  to  use  Eq.  (374)  as  a  model  for  the  numerical  algorithm.  A  fast  method  to  generate 
clouds  is  due  to  Perlin  [?].  We  instead  use  a  direct  implementation  since  we  are  not  limited  by  time 
considerations.  One  generates  a  cloud  s(x,y)  by  filtering  white  noise,  N(x,y),  with  a  filter  whose 
shape  in  the  frequency  domain  has  the  form  given  by  Eq.  (374), 


H' y  (Wj,  iOy  ) 


1 


(375) 


where  uiy)  is  the  transfer  function  of  the  linear  filter.  In  particular,  we  start  with  white  noise 

where  the  power  spectrum  is  constant 


Pn(^x,i^v)  =  Ao 


(376) 
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and  form  the  power  spectrum  of  the  cloud  by  way  of 


Psi^xi^y)  —  ^y)\  Pn(^xi  ^y)  (377) 

Then,  the  cloud  is  given  by 

s{x,  y)  =  j  j  hn(x  -x',y-  y')N(x',  y')dx'dy  (378) 

where  h^(x,y)  is  the  impulse  response  of  the  filter.  In  this  way  the  output  power  spectrum  of  the 
cloud  will  be  the  same  type  as  in  Eq.  (374). 

In  Fig.  1  we  show  several  clouds  generated  with  this  method.  In  every  of  the  five  pictures  we 
change  the  value  of  the  spectral  amplitude  7,  and  we  obtain  different  “types”  of  clouds.  In  particular 
in  Fig.  la  we  have  7  =  0,  that  is  we  generate  white  noise  without  any  filtering.  Then  in  Fig.  lb  with 
7  =  0.5  we  begin  observing  a  less  uniform  pattern  in  the  generated  image,  and  in  Fig.  lc  and  Fig.  Id 
obtained  respectively  with  7  =  1  and  7  =  1.5  we  have  the  more  cloud-like  aspect  of  natural  clouds. 
Fig.  le  with  7  =  2  gives  an  image  resembling  natural  clouds  quite  well.  Fig.  If,  with  7  =  5,  shows 
that  at  the  limit  the  generated  images  are  very  clustered. 
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Fig.  le.  7  =  2. 


Fig.  If.  7  =  5. 


30.2  Relative  phase-amplitude  importance 

We  write  the  Fourier  transform  in  terms  of  its  amplitude  and  phase 

F(ux,Uy)  =  A(ux,  (379) 

To  study  the  relative  importance  of  spectral  amplitude  and  phase  we  use  the  idea  of  Oppenheim  and 
Lim  and  reconstruct  the  image  using  phase  and/or  amplitude  only  to  see  the  relative  importance. 
In  fact  we  take  a  somewhat  more  general  approach.  We  define  a  modified  Fourier  transform  by  the 
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following 


(380) 


Fap{ujx,ujy)  =  Aa{ux,uy)e^x^ 
and  reconstruct  the  modified  image  by 

sap{x,y)  =  J  j  Fap(u)x,Uy)elx“x+iyUydu>xd(jjy  (381) 

In  Equation  (381)  a  and  (3  are  parameters  and  by  varying  a  and  /3  we  can  ascertain  the  relative 
importance  of  phase  and  amplitude.  In  particular  if  we  take  a  =  0  and  (1=1  then  we  are  reconstructing 
phase  only  images,  and  if  we  take  (3  =  0  and  a  =  1  then  it  is  amplitude  only  images.  Our  aim  is 
to  take  a  particular  cloud  and  reconstruct  it  by  the  above  procedure  with  a  variety  of  values  of  a 
and  (3  with  the  expectation  that  visual  inspection  will  give  us  some  indication  where  the  cloud-like 
information  resides,  amplitude  or  phase. 

In  Fig.  2  we  show  amplitude  only  and  phase  only  reconstruction  of  the  first  five  clouds  of  Fig. 
1.  In  particular,  in  the  first  column  of  the  figure  are  the  original  images.  Column  two  represents 
amplitude  only  reconstruction  and  column  three  represents  phase  only  reconstruction.  In  terms  of  the 
model  given  in  Sect.  30.1.  when  we  reconstruct  the  image  from  the  amplitude  we  take  a  =  1  and 
(3  =  0,  while  the  image  reconstructed  from  the  phase  content  has  a  =  0  and  (3=1.  From  a  study  of 
Fig.  2,  it  is  reasonable  to  conclude  that  both  amplitude  and  phase  are  important  for  the  reproduction 
of  the  original  image. 

In  Fig.  3  we  have  taken  the  image  of  the  realistic  cloud  case,  namely  the  one  where  7  =  2  in 
Fig.  1,  and  varied  a  and  (3  as  described  in  Section  3.  In  the  first  column  the  original  phase  is  kept 
(/3  =  1)  and  the  amplitude  significance  is  varied  according  to  a  =  0,  0.2,  0.4,  0.6,  0.8,  0.95.  In  the 
second  column  the  amplitude  is  kept  fixed  (a  =  1)  and  the  relative  importance  of  the  phase  is  varied 
by  taking  (3  =  0,  0.2,  0.4,  0.6,  0.8,  0.95.  Also  Fig.  3  proves  that  both  amplitude  and  phase  are 
important  to  generate  realistic  clouds. 
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Fig.  2 


a  =  0,  (3  =  1 
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30.3  The  1/f  model  for  cloud  generation 

In  the  previous  sections  we  have  considered  a  cloud  as  a  signal  and  we  have  investigated  where  in  the 
Fourier  transform  is  the  cloud-like  information  stored,  that  is  in  the  spectral  amplitude  or  spectral 
phase.  We  will  now  study  the  generation  of  clouds  from  a  different  perspective,  that  is  we  generalize 
the  algorithm  that  we  have  used  before  by  using  a  stochastic  fractional  differential  equation  model. 
In  particular  we  study  these  issues  on  artificially  generated  clouds. 
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It  is  generally  argued  that  clouds  have  an  1//  spectra,  and  it  has  been  well  known  that  this  can 
be  used  to  generate  clouds.  We  will  investigate  whether  this  is  indeed  true,  in  the  sense  that  do  other 
spectra  also  generate  clouds,  or  must  there  be  an  1//  behavior  in  the  spectra? 

The  1  //  approach  to  generate  clouds  assumes  that  the  power  spectra  of  the  generated  cloud  is  of 
the  form 


P(^xi  ^y)  — 


,27 


(382) 


[yjul  +  uy) 

where  7  is  the  1//  parameter  that  controls  the  cloud  like  appearance  of  the  image.  We  generated  a 
cloud  s(x,y )  by  filtering  white  noise,  N(x,y),  with  a  filter  whose  shape  in  the  frequency  domain  has 
the  form 


H'y  ) 


1 


1 

u/y 


where  H^((jOx,iOy)  is  the  transfer  function  of  the  linear  filter  and 


(383) 

(384) 


ui  = 


(385) 


In  particular,  we  started  with  white  noise  where 


PN  (Wj,,  iOy  )  -  N() 


(386) 


and  formed  the  power  spectrum  of  the  cloud  by 


Psi^xi^y)  —  ^y)  |  iVK,Wy) 

=  N0 

{s/ul+Vy)  7 

The  cloud  is  then  given  by 


(387) 

(388) 


s(x,  y)  =  JJ  hn(x  -  x',  y  -  y')N(x',  y')dx'dy'  (389) 

where  h^(x,y)  is  the  impulse  response  of  the  filter.  Eq.  (387)  has  the  power  spectrum  given  by  Eq. 
(374). 

In  Fig.  4  we  reproduce  some  clouds  generated  by  this  method  with  various  parameters.  We  point 
out  that  for  every  picture  we  are  using  a  two-dimensional  transfer  function  that  has  a  radial  symmetry, 
and  H{lo)  given  in  Eq.  (390)  is  precisely  the  value  of  that  transfer  function  at  any  angle  from  the 
reference  point  of  the  image  (located  in  the  middle  of  the  image),  and  at  the  radius  value  to. 
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30.4  Other  Power  Spectra  for  the  Generation  of  Clouds 


Our  aim  is  to  investigate  whether  indeed  only  1//  power  spectra  generates  cloud  like  images.  Towards 
that  end  we  first  consider  spectra  formed  by  the  process  described  above  but  where  we  take  the  transfer 
function  to  be  of  the  form 


H  M 


1 

(ico)' i'  +  ( iuj)s 


(390) 


This  produces  a  power  spectrum  given  by 


Psi^xi  bJy)  — 


(icvp  +  ( ilo)s 


1 

~~  U2^  +  C02S  +  [(-l)5  +  (-1)T')0'+'5wT'+<5 

1 

a;2')'  +  uj2S  +  2uH+l5  cos  [vr/2(7  —  5)] 
We  also  consider  transfer  functions  of  the  following  form 


(391) 

(392) 

(393) 


H{w) 


1 

uH  +  ico5 


in  which  case  the  power  spectrum  is 

P.s  (^ x  i^y)  — 

(J2'!'  +  <jJ25 


1 

uf<  +  ico5 

1 


(394) 


(395) 

(396) 


We  emphasize  that  this  power  spectrum  is  not  the  power  spectrum  of  the  sum  of  two  1/f  power 
spectra.  If  7  and  6  differ  by  multiples  of  ir,  then  clearly  Eq.  (391)  and  Eq.  (395)  are  the  same  power 
spectrum.  To  generate  the  cloud  we  use  the  procedure  described  in  Sect.  30.5.  Namely,  we  generate 
white  noise  and  construct  the  cloud  as  before,  but  we  use  for  the  transfer  function  Eq.  (390)  and  Eq. 
(394).  In  the  next  section  we  give  the  results  we  obtained  by  varying  7  and  6. 


30.5  Differential  Equation  Approach 

We  point  out  that  one  can  formulate  these  types  of  processes  and  procedures  as  fractional  differential 
equations.  Consider  the  following  fractional  differential  equation 

L>7x(f)  +  Dsx{t)  =  F(t)  (397) 


where  D 7  is  called  the  fractional  derivative  when  7  is  not  an  integer.  There  are  many  possible 
definitions  of  fractional  derivative.  Here  we  use  the  Weil  definition  [?] 

rt 


D^x{t)  = 


1 


r(7)  j- c 


(t  —  t'y  1x{t')dt' 


(398) 
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When  7  is  an  integer,  say  7  =  n,  this  reduces  to  the  n-the  derivative  of  x(t).  In  Eq.  (397)  F(t)  can 
be  any  function  but  here  we  take  it  to  be  white  Gaussian  noise. 

Now  consider  the  Fourier  transform  of  both  sides  of  the  differential  equation.  We  define  the 
Fourier  transform  G(uj)  of  a  function  g{t)  by 


GV)  =  F\g{t)\  =  -j=  J  g{t)t 

Also  we  define  X{oj)  to  be  the  Fourier  transform  of  x(t) 


e~iutdt 


X(u)=F[x(t)} 


Now  one  can  show  that 


F[D7x(t)]  =  (iu))1  F  [x(t)\  =  (iu)1  X(tjj) 
and  therefore  taking  the  Fourier  transform  of  both  sides  we  have  that 

F[D^x(t)  +  Dsx(t)}  =  F[F(t )] 


or 


giving 


(iuyx(  u)  +  (■ iuj)sX{u )  =  F[F{t) } 


X(u)  =  - 


F[F(t)} 


(iio)1  +  (iio)5 
Thus  we  see  that  the  transfer  function  is  equal  to 

H{u)  =  -  1 


(399) 

(400) 

(401) 

(402) 

(403) 

(404) 

(405) 


(iuj)' 1  +  ( iu)5 

that  is  Eq.  (390). 

For  the  transfer  function  given  by  Eq.  (394)  consider  the  following  fractional  differential  equations, 


i_7D7x(i)  +  =  F(t) 


For  this  equation  using  the  same  procedure  as  above  we  have  that 


xiu)  =  M 

v  u'y  +  ii05 


(406) 


(407) 


We  note  that  for  fractional  derivatives,  just  as  for  ordinary  derivatives,  we  have  that 


D 7  [ax(t)]  =  aD1x{t)  (408) 

where  a  is  an  arbitrary  constant. 

We  first  discuss  the  results  for  the  first  transfer  function,  Eq.  (390),  which  produces  the  power 
spectra  given  by  Eq.  (391).  In  Figs.  2a,  2b  and  2c  we  show  the  pictures  obtained  by  using  all  the 
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combinations  of  7  and  5  ranging  in  the  values  0.01,  0.5,  1,  1.5  and  2.  The  symmetric  combinations 
like  (7, 6)  =  (1,  2)  and  (7, S)  =  (2, 1)  have  been  avoided,  since  they  clearly  generate  identical  fractional 
equations  and  hence  processes  x(t),  as  can  be  seen  from  Eq.  (397).  The  rows  show  the  values 
7  =  2, 1.5, 1,  0.5, 0.01  (in  descending  order),  while  the  columns  are  associated  to  6  =  2, 1.5, 1,  0.5,  0.01 
(from  left  to  right).  Similar  results  have  been  obtained  in  Figs.  3a-3d,  where  we  have  used  the  transfer 
function  given  by  Eq.  (394).  All  the  possible  combinations  of  7  and  5  used  in  the  previous  case  have 
been  considered,  except  7  =  6  =  0.01,  which  corresponds  to  white  noise,  already  shown  in  Fig.  2c. 
Clearly  there  are  many  values  of  the  parameters  that  generate  cloud  like  pictures  and  hence  we  have 
to  conclude  that  it  is  not  the  case  that  clouds  must  have  an  1//  spectra. 

We  now  address  the  general  issue  of  1//  spectra.  It  has  often  been  said  that  many  things  in  nature 
have  1//  spectra,  examples  being  clouds,  speech,  natural  scenes,  etc.  and  that  this  is  a  fundamental 
characteristic.  The  definition  of  1//  spectra  often  gets  extended  to  mean  l//7.  We  are  currently 
studying  the  issue  as  to  what  is  common  about  the  spectra  that  produce  cloud  like  images.  We  specu¬ 
late  that  what  is  going  on  is  the  following.  First,  we  point  out  that  if  the  autocorrelation  function  is  of 
such  a  nature  that  there  is  no  correlations  at  all  we  get  white  noise  and  if  the  autocorrelation  function 
has  infinite  correlation  we  get  a  uniform  gray.  We  believe  that  cloud  like  images  are  formed  when  the 
autocorrelation  has  a  certain  range  and/or  spectral  shape  and  further  that  for  such  autocorrelation 
functions  there  are  many  types  of  spectra  that  can  produce  them,  not  only  1//.  Of  course,  there  is 
a  one  to  one  relationship  between  a  given  autocorrelation  function  and  power  spectra  but  what  we 
are  saying  is  that  for  a  class  or  range  of  correlation  functions  there  are  many  power  spectra  of  very 
different  functional  forms  that  can  produce  such  autocorrelation  functions. 

Insert  pictures  from  file  CloudFigures.tex  here 


30.6  Nonstationary  clouds 


Due  to  the  different  local  conditions  of  wind,  humidity,  and  temperature,  there  is  an  evident  variation 
in  the  uniformity  of  the  shape  of  clouds.  Since  clouds  vary  in  nature  with  position,  it  is  natural  to 
think  that  also  their  spectral  representation  will  be  a  function  of  position.  It  is  hence  interesting  to 
build  a  realistic  model  of  such  nonstationarities. 

We  have  seen  that  the  power  spectrum  Ps(LOx,u>y)  of  a  cloud  s(x,y )  can  be  modeled  as 


P.s  {f^xi  W  y  ) 


1 


(409) 


To  obtain  this  power  spectrum  we  have  processed  white  Gaussian  noise  f(x,y )  with  a  filter  that  has 
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a  transfer  function  given  by 


H  ( ivx  ,  /)  — 


/wg  +  wj* 


We  remind  that  if  we  start  with  white  Gaussian  noise,  which  has  a  constant  power  spectrum  given  by 


the  result  of  the  filtering  is  [?] 


Pfi^xi^y)  —  N0 


Psif^xi^y)  —  \P{^xi^y)\  Pfi^xi^y) 
=  No 

(M+<4) 7 


This  is  exactly  the  type  of  spectrum  we  seek,  namely  Eq.  (409).  In  the  spatial  domain  the  filtering  is 
given  by  a  two-dimensional  convolution 


s(x,  y )  =  j  J h(x  -x',y-  y')f{x ',  y')dx'dy' 


where  h(x,  y )  is  the  impulse  response  of  the  filter,  obtained  from  the  transfer  function  via  Fourier 


inversion 


h(x,y )  =  JJ  H(ux,uy)e™xX+tWyVdujxdujy 


The  type  of  clouds  generated  by  the  described  method  depends  on  the  parameter  7.  In  Fig.  30  we 
show  the  cloud  corresponding  to  7  =  1,  while  in  Fig.  31  the  case  7  =  1.5  is  shown.  One  notices  that 
as  7  increases,  the  clouds  are  bigger  and  change  slower  with  respect  to  position.  This  effect  is  due  to 
the  increased  lowpass  nature  of  the  filter  H(<jjx,u)y)  used  to  generate  the  cloud. 

Since  the  clouds  vary  in  nature  with  position,  it  is  natural  to  think  that  also  their  spectral  repre¬ 
sentation  will  be  a  function  of  position.  Given  that  we  are  using  a  1  //7  spectral  model,  an  immediate 
way  of  modeling  a  nonstationary  cloud  is  to  make  7  a  function  of  position,  that  is  to  have  7 (x,y). 
In  order  to  do  this,  we  use  the  filter  approach  reviewed  in  the  introduction  where  now  we  take  the 
transfer  function  to  depend  on  position  in  the  following  way 


H(ux,uy,x,y)  = 


(v 


That  is,  the  dependence  of  the  filter  shape  on  position  is  obtained  by  making  7(2;,  y)  a  function  of 
position.  Therefore  the  convolution  operation  of  Eq.  (414)  that  generates  the  cloud  s(x,  y)  from  the 
white  Gaussian  noise  f(x,y )  will  be  replaced  by  the  general  integral 

s{x,y)=  [  fh(x,x,,y,y')f(x',y')dxldy'  (417) 
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Figure  30:  Stationary  cloud  obtained  with  7=1. 


Figure  31:  Stationary  cloud  obtained  with  7  =  1.5. 
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(418) 


where  now  h(x,  x' ,  y,  y')  is  the  Green’s  function  given  by 

h{x,x',y,y')  =  H(ux,uy,x,y)el“xX'+lUyy'  duxduy 

The  corresponding  algorithm  is  a  direct  implementation  of  the  Green’s  integral  of  Eq.  (417). 

We  now  test  the  method  by  imposing  different  laws  for  the  parameter  7 (x,y). 

Radial  case 

We  first  impose  7  to  be  constant  for  all  the  positions  that  have  the  same  distance  from  a  given 
reference  point.  We  choose  as  the  reference  point  the  position  (x,  y)  =  (0,  0),  and  we  model  7  with 

7(r)  =  7l  +  72  R  71  r  (419) 

where  r  is  the  radius  given  by 

r  =  \/x2  +  y2  (420) 

and  R  is  the  maximum  distance  from  the  reference  point,  that  is  reached  in  the  corner  {xmaXiVmax) 
that  is  the  farthest  from  the  reference  point 

R  =  \JX'MAX  +  y'l/IAX  (421) 

Hence  Eq.  (419)  imposes  a  linear  variation  of  7  from  an  initial  value  of  71  in  the  reference  point  to  a 
final  value  of  72  in  the  farthest  point.  Fig.  32  represents  the  values  taken  by  7  on  the  image,  when 
71  =  1  and  72  =  1.8.  In  Fig.  33  we  show  the  result  of  our  method.  We  notice  that  the  uniformity 
of  the  clouds  vary  with  position.  In  the  lower  bottom  corner  the  clouds  are  more  similar  to  the  ones 
of  Fig.  30,  obtained  for  the  stationary  case  with  7  =  1.  While  we  reach  the  upper  top  corner  we  see 
that  the  clouds  become  similar  to  the  one  of  Fig.  31,  that  represents  the  stationary  case  with  7  =  1.5. 
We  have  hence  obtained  our  goal  to  generate  nonstationary  clouds  by  changing  the  local  shape  of  the 
filter  H(ux,uy,x,y). 

Lateral  case 

We  now  impose  7  to  be  a  function  of  x  only,  that  is 

7 (x)  =  71  +  72  ~7lx  (422) 

xmax 

Hence  we  have  chosen  7  to  increase  linearly  from  71  to  72  with  the  x  position,  and  to  be  constant  with 
respect  to  y.  Fig.  34  shows  the  values  taken  by  7,  when  71  =  1  and  72  =  1.8.  In  Fig.  35  the  clouds 
obtained  in  this  case  are  shown.  We  notice  that  the  clouds  located  towards  the  left  side  of  the  picture 
are  similar  to  the  ones  obtained  in  the  stationary  case  with  7  =  1,  shown  in  Fig.  30.  While  we  move 
towards  the  right  side  of  the  picture  we  see  that  the  clouds  become  more  similar  to  the  ones  obtained 
in  the  stationary  case  with  7  =  1.5  and  represented  in  Fig.  31. 
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Figure  32:  These  are  the  values  taken  by  the  7  parameter  in  the  model  of  Eq.  (  418).  The  parameter 
7  follows  Eq.  (  419),  with  71  =  1  and  72  =  1.8. 
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Figure  33:  This  picture  shows  nonstationary  clouds  that  have  been  obtained  by  imposing  the  parameter 
7  to  change  as  per  Eq.  (419).  The  clouds  in  the  region  around  the  bottom  left  corner  are  similar  to 
the  ones  of  Fig.  30,  while  the  clouds  around  the  upper  right  corner  are  better  represented  by  the  ones 
in  Fig.  31. 
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Figure  34:  These  are  the  values  taken  by  the  7  parameter  in  the  model  of  Eq.  (418).  The  parameter 
7  follows  Eq.  (  422),  with  71  =  1  and  72  =  1.8. 
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Figure  35:  This  picture  shows  nonstationary  clouds  that  have  been  obtained  by  imposing  the  parameter 
7  to  change  as  per  Eq.  (422).  The  clouds  in  the  region  around  the  left  side  of  the  image  are  similar 
to  the  ones  of  Fig.  30,  while  the  clouds  in  the  region  around  the  right  side  of  the  image  are  better 
represented  by  the  ones  in  Fig.  31. 
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31  Adjustable  bandwidth  concept 


The  adjustable  bandwidth  concept  (ABC)  is  a  method  that  enhances  the  performances  of  detection 
and  estimation  algorithms  [?,  ?,  ?].  The  aim  is  to  enhance  nonstationary  signals  in  noise,  that  is, 
to  bring  out  the  main  features  of  signals  that  may  be  buried  in  noise.  This  is  essential  for  detection 
and  classification  and  the  method  may  be  thought  of  as  a  “preprocessing  algorithm”  for  classification 
applications  and  in  particular  for  automatic  target  recognition  methods.  In  this  article  we  show  that 
the  ABC  approach  is  a  time-frequency  method  in  that  its  basic  idea  is  to  treat  time  and  frequency 
jointly.  In  particular  we  show  that  the  algorithm  transforms  a  time-frequency  distribution  by  means  of 
individual  time- frequency  transformations  into  a  set  of  final  time- frequency  distributions.  Each  step 
of  the  procedure  can  be  understood  and  formulated  in  terms  of  the  kernel  method,  a  technique  that 
is  standard  in  time-frequency  analysis  [?,  ?].  Once  the  procedure  is  formulated  in  terms  of  kernels, 
the  ABC  properties  can  be  better  studied  and  can  suggest  further  enhancement  and  effectiveness  of 
the  method.  We  develop  an  exactly  solvable  example  that  illustrates  some  of  the  main  issues  and  also 
shows  the  behavior  of  the  method. 

31.1  Time-frequency  and  the  kernel  method 

A  fundamental  idea  of  time-frequency  distributions  is  the  concept  of  a  kernel.  The  kernel  character¬ 
izes  the  distribution  and  its  properties  and  typically  by  examining  the  kernel  one  can  ascertain  the 
properties  of  the  distribution.  There  have  been  many  methods  and  distributions  proposed  over  the 
years,  among  them  the  Wigner  distribution,  the  spectrogram,  the  Choi- Williams,  Margenau-Hill  or 
Rihaczek  and  the  Zam  distributions  and  several  others.  Among  the  many  areas  to  which  they  have 
been  applied  are  biomedical  signal  analysis  (e.g.,  heart  sounds,  heart  rate,  the  electroencephalogram 
(EEG),  the  electromyogram  (EMG)  and  others),  machine  fault  monitoring,  radar  and  sonar  signals, 
acoustic  scattering,  wave  propagation,  speech  processing,  analysis  of  marine  mammal  sounds,  musical 
instruments,  linear  and  nonlinear  dynamical  systems,  among  many  others.  The  primary  reason  for 
the  applicability  of  these  methods  to  such  a  variety  of  fields  is  that  in  all  the  cases  cited,  the  spectra 
of  the  signals  of  interest  change  with  time  and  these  changes  are  fundamental  to  understand  as  they 
reflect  the  source  and/or  propagation  medium  [?]. 

All  time- frequency  distributions  may  be  generated  from  the  general  class  that  is  given  by  [?] 

C \t,u)  =  ^  fff  s*(u  ~  \T^  s(«+  lT)  <Kd,T)  e-jdt-jTUJ+jdududrde  (423) 

where  c t>{9 ,  r)  is  a  two  dimensional  function,  the  kernel.  An  alternative  and  useful  formulation  is  to 
write 

C(t,w)  =  4\2  ff  M(e,T)e-jet~jTUJd6dT  (424) 
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where 


M(0,t )  =4>(6,t)  j  s*(u  —  -t)  s(u  +  -t)  e^9udu 
=  (/)(0,t)A(0,t) 


and  where  A(0,  r)  is  the  standard  ambiguity  function  commonly  used  in  radar  for  the  design  of  signals. 
The  function  M(6,t)  is  called  the  characteristic  function  of  the  distribution. 

We  now  briefly  mention  some  kernels  and  their  respective  distributions.  If  the  kernel  is  taken 
to  be  one  then  the  Wigner  distribution  is  obtained  [?].  To  avoid  the  cross  term  problems  of  the 
Wigner  distribution  a  kernel  methodology  has  been  developed  for  kernel  design  and  these  ideas  were 
initiated  by  Choi  and  Williams  and  Zhao,  Atlas  and  Marks  [?].  The  Choi-Williams  kernel  is  given 
by  <!>(6,t)  =  exp (— #2t2/ct)  and  the  Zam  kernel  by  cj)(6,  r)  =  g(r)  |r|  .  Williams  and  co-workers 

devised  and  crystallized  the  idea  of  kernel  design  [?,  ?].  They  developed  a  methodology  for  the 
construction  of  distributions  with  desirable  properties.  Of  particular  interest  is  the  spectrogram, 
given  by 


Cs(t,u)  = 


w(t  —  t')s{t')e  lut  dt' 


where  w(t)  is  the  window.  The  kernel  for  the  spectrogram  is  given  by 

<j>(0,T)  =  J  w* (u  —  -t)  iu(u  +  -t)  e^9udu  (428) 

That  is,  the  kernel  of  the  spectrogram  is  the  ambiguity  function  of  the  window.  We  point  out  that 
the  squared  magnitude  of  the  wavelet  transform  can  be  used  as  a  time  frequency  distribution  and 
formulated  in  the  above  terms  [?,  ?].  Also,  a  number  of  fundamental  results  on  the  relation  between 
the  distribution  and  kernel  have  been  developed,  over  the  last  ten  years,  the  paper  by  Loughlin  et.  all 
being  fundamental  [?] 

While  the  kernel  method  has  generally  been  used  to  study,  obtain,  and  characterize  distributions, 
we  discuss  here  the  fact  that  it  can  also  be  used  as  a  way  of  transforming  distributions.  The  reason  we 
emphasize  this  is  that  the  ABC  method  is  a  sequence  of  time-frequency  distribution  transformations. 
Here  we  give  some  general  properties  of  transformation.  Suppose  we  have  two  distributions,  C\  and 
C2,  with  corresponding  kernels,  <f>i  and  <f>2.  Their  respective  characteristic  functions  are  related  by 

Mi(6,t)  =  (429) 

<h{0,T) 

This  relationship  connects  the  characteristic  functions  of  any  two  distributions.  To  obtain  the  rela¬ 
tionship  between  the  distributions  one  uses  Eq.  (424)  to  obtain 

Ci(t,u)  =  [[  gi2(t'  —  —  u)  C2{t' ,  Lu')dt'd(jj'  (430) 
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with 


*““■") Vu+>T“d0dT  <431) 

What  this  shows  is  that  one  distribution  can  be  transformed  into  another  by  way  of  a  two-dimensional 
convolution  where  the  convolution  function  is  given  in  terms  of  the  kernels  of  each  distribution. 

31.2  ABC  Method 

The  ABC  approach  is  a  multi-stage  procedure  that  operates  on  a  starting  time- frequency  distribution, 
Cs(t,u)),  and  generates  N  output  distributions,  C\(t,  lo),  C2{t,u),  ...  ,  Cjv(t,u;).  In  this  paper  we 
show  that  the  procedure  at  each  stage  is  a  time-frequency  distribution  transformation  and  that  the 
end  distribution  can  be  characterized  by  a  kernel  which  is  a  functional  of  the  kernels  at  each  stage. 

We  first  give  an  overall  view  of  the  basic  algorithm  and  indeed  we  present  two  but  equivalent 
views.  The  first  formulates  everything  in  the  time- frequency  plane  and  the  second  formulates  some 
of  the  steps  in  the  kernel  or  ambiguity  function  domain  as  discussed  above.  Each  formulation  has 
certain  advantages,  the  first  formulation  allows  a  clearer  mathematical  view  of  the  issues  but  the 
second  formulation  lends  itself  to  more  detailed  discussion  of  the  implementation  issues.  We  present 
both  formulations  but  in  this  paper  we  discuss  in  detail  the  second  one.  The  first  view  is  presented  in 
detail  in  [?]  where  also  an  analytic  example  is  given.  The  equivalence  of  the  two  formulations  is  given 
in  [?]. 

In  the  first  formation  the  basic  idea  is  to  take  a  time-frequency  distribution  and  from  it  obtain  N 
new  distributions.  The  new  distributions  have  progressively  different  perspectives  or  resolutions  of  the 
original  distribution.  For  the  starting  distribution  we  use  Cs(t,u ;),  and  the  N  resulting  distributions 
are  C2(t,co),  ...  CN(t,u).  However  the  actual  iterative  process  is  done  on  N  intermediate 

distributions  that  we  call  ^(f,  w), ...,  Ajy(t,cu).  The  iterative  process  is 

Am(t,u)  =  Am-i(t,u)  -  J  hm-i(t  -  t',u  -  J)  Am-i(t' ,u')dt' dJ  (432) 

and  to  start  one  takes  Ai(t,u)  =  Cs(t,uj).  The  hm(t,ui )  are  fixed  functions.  We  note  that  at  each 
stage  one  removes  from  the  distribution  the  effect  that  the  second  term  of  Eq.  (432)  produces.  The 
iterative  scheme,  Eq.  (432),  is  self  contained  and  produces  distributions  with  certain  properties  that 
are  achieved  by  the  judicious  choice  of  the  hm(t,iv)  so  chosen  to  achieve  required  results,  such  as 
frequency  smoothing.  Now,  to  further  manipulate  and  enhance  the  result  we  obtain  Cm(t,  oj)  from 
Am(t,u)  by 

Cm(t,u)  =  J  gm(t  -  t',u  -  J)Am(t' ,u')dt'du'  (433) 

where  again  gi(t,  oj),  <72 (i,  w), gjv(i,  to)  are  fixed  filter  functions  prechosen  to  achieve  desired  char¬ 
acteristics.  Of  course  one  can  just  stick  with  the  A’s  and  from  a  mathematical  point  of  view 
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that  can  be  achieved  by  taking  delta  functions  for  the  g’s,  in  which  case  the  result  would  be  that 
Cm(t,uj)  =  Am(t,uj).  Also,  one  can  combine  Eq.  (433)  and  Eq.  (432)  but  we  have  found  it  clearer 
to  write  the  procedure  in  the  above  form.  We  note  that  both  transformations  in  Eq.  (433)  and  Eq. 
(432)  are  of  the  form  as  given  by  Eq.  (430). 

In  the  second  formulation,  at  each  stage  of  the  algorithm  there  is  an  input  distribution  Am(t,Lu), 
an  intermediate  distribution  and  an  output  distribution  Cm(t,tv).  The  input  distribution  at 

the  stage  1  is  A\(t,u)  =  Cs(t,iv).  The  distributions  Am(t,uj)  and  Bm(t,Lo)  are  needed  to  calculate 
the  output  and  to  go  on  to  the  next  stage.  In  particular  we  tabulate  here  the  overall  steps: 

Am(t,uj)  =  Am_i(t,w)  —  Bm—i(t,Lu)  (434) 

Bm{t,uj )  =  is  calculated  from  Am(t,Lo)  through  a  frequency  average  (435) 

Cm(t,u ;)  =  is  calculated  from  Bm(t,  w)  through  a  time  average  (436) 

In  addition,  we  will  need  the  characteristic  functions  for  each  distribution  and  we  use  the  following 

notation  am(9,  r),  bm(9,r)  and  cm(0,  r).  That  is 

am(0,  r)  =  jj  Am(t,  uj)  e*et+iTbJdtdujT  (437) 

and  similarly  for  bm(0,T),  and  cm(0,r).  We  now  describe  the  two  steps  indicated  above.  For  both  of 
them  we  will  use  the  characteristic  function  formulation,  because  it  makes  the  derivation  much  easier. 

31.3  First  sub  step:  Calculation  of  Bm(t,cu )  from  Am(t,u)  through  a  frequency 
average 

The  distribution  Bm(t,uj )  is  obtained  from  Am(t,Lo)  through  a  lowpass  filtering,  with  the  aim  of 
smoothing  out  noise  along  the  frequency  axis  of  the  input  distribution  of  the  stage.  The  amount  of 
filtering  is  in  general  a  function  of  the  stage.  In  the  characteristic  function  domain  we  have  that,  at 
stage  m  =  1,2 , ...  IV 

bm(9,  t)  =  V2tt Hm(T)am(e,T)  (438) 

where  Hm{r)  is  the  transfer  function  of  the  lowpass  filter  at  the  m-th  stage.  We  have  defined 

Hn(t)  =  1  (439) 

Also,  because  of  Eq.  (434)  we  see  that  for  m  =  2,  3, . . .  N 

am{0,  r)  =  [1  -  V2nHm-i(T)]am-i(8 ,  r)  (440) 

m— 1 

=  II  [1  -  V2 ^m_Kr)]Ms(0,  r)  (441) 

l=i 
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Putting  together  Eq.  (438)  and  Eq.  (440)  we  obtain 


bm(6,r)  =  </w(t)Ms(0,t) 


m—  1 


\/2i xHm(r)  [1  -  V2tt Hm_i(r)\ 


i=i 

where  we  have  defined  the  frequency  averaging  kernel  as 

m— 1 


Ms(9,t) 


(frujmi'r)  =  \/2i xHm(T)  [1  -  v/2vrpfm_i(r)] 


i=i 


(442) 

(443) 


(444) 


31.4  Second  sub  step:  Calculation  of  Cm(t,Lv)  from  Bm(t,ou )  through  a  time  average 

The  distribution  Cm{t,io)  is  obtained  from  Bm(t,  uj)  through  a  moving  average  in  time,  with  the  aim 
of  reducing  noise  along  the  time  direction  of  the  input  distribution  of  the  stage.  If  we  consider  the 
m-th  stage  of  the  algorithm,  m  =  1,  2, . . . ,  N,  the  time  averaging  can  be  written  in  the  time-frequency 
domain  as 

Y  rt+Tm/ 2 

Cm(t,uj)  =  —  /  Bm(t' ,  iv)dt'  (445) 

^  in  Jt-Tm/ 2 

=  [ PTm{t-t')Bm{t' ,u)dt'  (446) 

p  m  J 

=  ^~Prm(t)*Bm(t,u)  (447) 

-L  m 

where  the  star  sign  indicates  convolution  and  PTm(t )  is  the  rectangular  window  defined  as 


PbM  = 


1 

0 


-Tm/2  <  t  <  Tm/2 
elsewhere 


(448) 


and  also  &at(0,t)  =  ajy(9,T).  In  the  characteristic  function  domain  we  have  that 


cm(9,r)  =  <j>tm{9)bm(9  ,t)  (449) 

=  \f2 n^PTrn{9)bm(9,T)  (450) 

where  PTm(9 )  is  the  inverse  Fourier  transform  of  PTm{t)  and  we  have  defined  the  time  averaging  kernel 
as 

<t>tm{9)  =  T^-PTrn{9)  (451) 

-*■  m 
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31.5  Time- Frequency  averaging  in  the  ABC  algorithm 

We  will  now  derive  the  exact  expression  for  the  output  distributions  generated  by  the  algorithm.  We 
see  by  considering  Eq.  (449)  and  Eq.  (442)  that  the  output  of  stage  m  is 


Cm{9 1  T ) 


^W«m(r)Ms(0,r) 


-  771—1 

27 T—PTrn(0)Hm(T)  [1  -  V2 nH,m_i(T)} 

J-  m  ,  n 


<f>m(0,T)Ms(O,T) 


(452) 

(453) 

(454) 

(455) 


where 


<M0>r)  =  &mW«m(r)  (456) 

-  771—1 

=  2tt— PTm(9)Hm(r)  JI  [1  -  (457) 

m  i=i 

is  the  kernel  of  the  distribution  at  stage  to,  that  is  hence  independent  from  the  input  distribution.  We 
note  that  at  stage  m  =  1  the  multiplication  operation  has  to  be  set  equal  to  1. 

Convergence  and  weighted  averages .6  We  now  address  the  issue  of  the  convergence  of  the  algorithm. 
We  have  found  from  experience  that  the  algorithm  always  converges  and  we  believe  that  some  insight 
into  this  issue  can  be  gained  by  the  following  considerations.  In  [?]  it  is  shown  for  the  discrete-time 
case  that  the  ABC  frequency  averaging  of  the  input  is  equivalent  to  a  filter-bank  decomposition. 
Therefore,  just  as  the  output  components  of  a  filter  bank  can  be  recombined  to  “converge”  to  the 
input,  the  outputs  identified  by  Eq.  (433)  can  be  recombined  to  “converge”  to  the  input.  Also,  in  [?] 
a  simple  analytic  example  is  presented  and  at  least  for  that  example  one  can  see  the  that  indeed  at 
each  step  one  gets  the  appropriate  convergence.  One  of  the  advantages  of  formulating  the  algorithm 
by  way  of  Eqs.  (433)  and  (432)  is  that  it  crystallized  the  procedures  and  allows  for  the  possibility  of 
a  mathematical  study  of  the  convergence.  This  issue  is  currently  being  studied. 

In  the  above  some  parameters  are  estimated  and  averaged  simply  but  of  course  one  could  use 
weighted  average.  Weighted  average  is  always  an  option  and  any  specific  choice  of  filter  parameters 
will  always  be  “optimal”  for  some  signal  scenarios  and  “sub-optimal”  for  others.  This  paper  focuses  on 
formulating  the  ABC  process  in  continuous-time  and  continuous-frequency,  from  the  original  discrete¬ 
time,  discrete- frequency  formulation.  This  new  formulation  then  lends  itself  to  study  this  issue.  Similar 
issues  regarding  weighted  averages  are  considered  in  [?]. 

6We  thank  the  referees  for  bringing  some  of  these  points  to  our  attention. 
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Relation  between  direct  time-frequency  formulation  and  the  ambiguity  function  domain.  We  repeat 
here  the  two  formulations  of  the  ABC  algorithm  for  convenience.  In  time-frequency  we  have  the  two 
equations 


Am(t,uj)  =  Am_i(t,uj)  (458) 

—  j  hm-i(t  —  t',u  —  u/)  ui')dt' dui' 

Cm(t,u)  =  J  gm{t  -  t',uj  -  J)Am(t' ,J)dt’dJ  (459) 

with  Ai(t,ui)  =  Cs(t,uj).  In  the  characteristic  function  domain  it  is 

cm{0 ,  t)  =  t)Ms(0 ,  r)  (460) 

=  ^J0)^m(r)Ms(0,T)  (461) 

=  2n^-Prm(e)Hm(T)]  x  (462) 

-Lm 

~m—  1 

n  [1  -  Ms(d,r)  (463) 

.  l=i 


To  point  out  the  connection  we  write  Eq.  (458)  and  Eq.  (459)  in  the  characteristic  function  domain 

am(0,T )  =  [l  -  hm(6,  t )]  r )  (464) 

=  Wk=i  t1  -  hk(0,T)\  ai(6,r)  (465) 

and 

cm{0,  t)  =  gm(6,  t )am  (466) 

=  9m(0,  r)nr=i  i1  ~  hk(6,  r)]  ai(6»,  r)  (467) 

where  hm(6,  r)  and  gm(0,  r)  are  the  filters  hm(t,  ui)  and  gm{t,uj)  in  the  characteristic  function  domain. 
Since 

ai(6,  t)  =  Ms{6,  t)  (468) 

we  see  that  if  we  choose 

t)  =  gm(0,  r)n^Li  [l  -  hk(e,  r)]  (469) 

we  have  the  equivalent  formulation  of  the  characteristic  function  domain.  By  further  relaxing  the 
constraint  of  the  time  kernel  <f>tm(0)  to  be  a  function  of  6  only  and  the  frequency  kernel  0ajm(r)  to  be 
a  function  of  r  only  we  establish  the  connection 


4>tm{0,T)  =  gm{6,T)  (470) 

(r)  =  nr=i  [!  -  r)]  (471) 

thus  showing  the  one  to  one  relationship  among  the  two  formulations. 
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31.6  Examples 

To  demonstrate  the  application  of  the  ABC  process  we  present  a  number  of  cases.  The  first  is  a 
signal  composed  of  the  sum  of  a  sinusoid,  a  chirp  (swept  frequency),  an  impulse,  and  additive  white 
Gaussian  noise.  The  term  “white”  refers  to  the  fact  that  the  noise  power  is  uniform  across  frequency. 
For  this  case,  the  noise  power  is  low,  with  the  signal-to-noise  power  ratio  (SNR)  set  to  70  dB  for  the 
sinusoid  as  measured  in  a  frequency  bin.  The  spectrogram  generated  is  based  on  a  1024  sample  FFT, 
with  a  rectangular  data  window,  no  data  overlap  between  time  segments  and  logarithmic  scaling.  The 
amplitudes  of  the  chirp  and  impulse  signal  components  are  set  to  achieve  similar  power  levels  per 
FFT  bin  as  the  sinusoid  component.  The  intent  is  to  demonstrate  separation  of  these  unique  signal 
components  into  the  various  output  stages  by  the  ABC  process. 

The  second  example  is  essentially  the  same  as  the  first,  with  the  exception  that  the  additive  noise 
component  is  increased  substantially  such  that  the  sinusoidal  SNR  is  now  10  dB  per  FFT  bin.  Thus 
the  sinusoid,  chirp  and  impulse  are  the  same  as  in  the  first  signal.  For  each  case,  the  signal  is  quantized 
to  16  bits  for  storage  to  file  prior  to  the  ABC  processing. 

For  both  signals,  the  ABC  parameters  are  set  to  the  same  values.  In  particular,  3  processing  stages 
are  used.  The  first  stage  lowpass  filter  is  a  finite  impulse  response  (FIR)  filter  of  129  coefficients,  each 
coefficient  equal  to  1/129.  Likewise,  the  stage  2  lowpass  filter  is  a  7  coefficient  FIR  filter,  each 
coefficient  equal  to  1/7.  Thus  both  filters  are  implemented  as  simple  moving  averages.  Note  that  filter 
design  will  affect  the  ABC  processing  performance.  While  the  parameters  chosen  are  appropriate  for 
the  signal  scenarios  presented,  these  parameters  have  not  been  optimized.  (The  effect  of  this  choice 
of  parameters  will  be  apparent.)  In  addition  to  frequency  averaging,  time  averaging  is  accomplished 
in  the  second  and  third  stages  of  the  implemented  ABC  process.  An  average  of  2  time  segments  is 
accomplished  in  the  second  stage,  and  an  average  of  10  time  segments  is  accomplished  in  the  third 
stage. 

Results  are  shown  in  Figs.  1-4  for  the  70  dB  SNR  case,  and  Figs.  5-8  for  the  10  dB  SNR  case. 
As  seen  in  Fig.  1,  the  sinusoid  is  at  bin  240  for  all  segments.  The  chirp  component  sweeps  up  in 
frequency  from  bin  0  to  511,  then  back  down  to  0.  The  impulse  occurs  at  time  segment  256.  Note  that 
in  progressing  from  stage  to  stage,  the  desired  components  are  much  more  distinguishable  from  each 
other,  and  from  the  background  additive  noise.  However,  note  also  that  the  energy  of  the  chirp  signal 
is  present  in  both  stage  1  and  stage  3,  in  addition  to  the  dominant  chirp  output  in  stage  2.  Likewise, 
the  sinusoidal  component  is  dominant  in  stage  3,  but  also  noticeable  in  stage  2.  This  observation 
relates  to  the  need  for  optimization  of  the  filter  parameters  of  the  ABC  process,  when  there  exists 
a-priori  information  regarding  the  input. 

In  Fig.  5,  note  the  adverse  affects  of  increasing  the  additive  noise  component.  The  sinusoid,  chirp 
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Figure  36:  Input  distribution  Cs(t,u)  when  SNR  =  70dB. 


Figure  37:  Output  distribution  C\  (t.  u)  at  stage  1  ( SNR  =  70dB). 
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Figure  38:  Output  distribution  C^i,  w)  at  stage  2  ( SNR  =  70dB). 


Figure  39:  Output  distribution  C^(t,Lo)  at  stage  3  ( SNR  =  70dB). 


90 


Figure  40:  Input  spectrogram  Cs(t,uj)  when  SNR  =  lOdB. 


Figure  41:  Output  distribution  C\(t,Lo)  at  stage  1  ( SNR  =  lOdB). 
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Figure  42:  Output  distribution  C^i,  w)  at  stage  2  ( SNR  =  lOdB). 


Figure  43:  Output  distribution  C^(t,Lo)  at  stage  3  ( SNR  =  lOdB). 
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and  impulse  are  difficult  to  discern.  Even  so,  the  benefits  of  the  ABC  process  are  easily  observed.  The 
desired  signal  components  are  separated  in  preparation  for  post-processing,  such  as  threshold-based 
detection. 

Example  2 

We  consider  a  specific  example  to  see  the  effect  of  the  procedure 
with  no  amplitude  modulation 

sit)  =  ejkt 2/2 

For  the  starting  distribution  we  take 

Cs(t,ca )  =  6(u  —  kt)  (473) 

which  is  the  Wigner  distribution  of  the  signal)?,  ?].  The  instantaneous  frequency  is  c Oi  =  kt  and 
therefore  the  distribution  is  totally  concentrated  along  the  instantaneous  frequency.  For  this  example 
we  calculate  C\(t,Lu),  A2(t,  w)  and  (72(t,u;). 

For  the  filters  we  take 

gi(t,  uj)  =  ^  exp(— out2  -  /3iw2)  (474) 

7 r 

g2(t,u)  =  v/a2/7rexp(-a2t2)5(w)  (475) 

hi(t,u>)  =  \/ /ir5(t)  exp(— /?iw2)  (476) 

We  have  taken  these  filters  because  they  are  directly  connected  in  a  simple  way  to  Noga’s  original 
formulation  of  the  algorithm.  In  particular,  the  filter  gi(t,  ca)  is  the  combination  of  the  lowpass  filtering 
and  the  time  averaging  performed  by  the  algorithm  at  the  first  stage.  The  lowpass  filter  has  impulse 
response  h\{t,(v)  and  the  first  time  average  is  exp(— aif2)5(ti;).  The  time  average  for  the  second  stage 
is  given  by  <72(t,u;).  All  these  filters  are  defined  in  the  time-frequency  plane,  while  in  the  original 
formulation  they  are  defined  on  time  and  frequency  separately. 

A  straightforward  calculation  gives  the  following  results 


For  the  signal  we  take  a  chirp 


(472) 


aiPi 

(ctq  +  P\k2) 


exp 


j.2  a  ,  ,2  ,  (ait  +  Pikuf 
—  Qt\t  —  Pl^  i 


UlPl 


7r(ai  +  Pi  k2) 


exp 


—a-iPi 


+  P\k2 
(a;  —  kt)2 
a\  +  Pik 2 


(477) 


(478) 


A2  =  5(lu  —  kt)  —  \J Pi/ir  exp  —  P\{u)  —  kt)2 


(479) 
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(480) 


C2(t,u>)  =  —  v/a2/7rexp(-a2(w  -  kt)2/k2) 


ct2Pi 


ir(a2  +  Pi  k2) 


exp 


-a2r  -  Piuz  + 


2  (a-2t  +  Pikup 


a2  +  P\k2 


=  -rT\/Q!2/vi'exp(-a2(a;  -  kt.)2/k2) 
\k\ 


(481) 


a2P\ 


7r  (a2  +  Pik2) 


exp 


—a2P\ 


( u}  —  kt )2 
a2  +  P±k2 


We  choose  the  following  parameters:  aq  =  1000,  a2  =  100,  Pi  =  0.01  and  k  =  2n  and  we  show 
the  quantities  Cs{t,u>),  Ci(t,Lo)  and  C2(t,ui).  In  Fig.  44  we  plot  the  input  distribution  Cs{t,cu),  that 
is  a  chirp  concentrated  along  the  instantaneous  frequency  LOi  =  pt.  In  Fig.  45  we  have  C\{t,uj),  the 
output  of  the  first  stage  of  the  algorithm.  We  notice  that  the  chirp  has  been  spread  in  both  time 
and  frequency.  This  effect  is  generated  by  the  filter  g\ (t,  cu)  that  performed  a  time  and  frequency 
smoothing  on  the  input  distribution.  The  importance  of  the  smoothing  is  related  to  noise  suppression, 
but  here  we  investigated  a  noise  free  deterministic  case  in  order  to  understand  the  basic  operation 
performed  by  the  method.  Finally  in  Fig.  46  we  show  the  second  output  of  the  algorithm,  C2(t,Lo). 
With  judicial  choices  of  filter  parameters,  we  can  adjust  the  amount  of  frequency  and  time  averaging 
in  each  stage.  In  the  same  fashion  as  evaluating  time  filters  with  an  input  impulse  and  observing 
the  output  response,  here  we  have  chosen  a  swept  impulse  in  the  time-frequency  plane  to  evaluate 
and  observe  the  output  results.  We  observe  that  the  first  stage  output  responds  only  slightly  to  the 
chosen  input  signal,  while  the  second  stage  responds  substantially.  This  is  due  to  the  fact  that  the 
input  signal  is  narrow  in  (instantaneous)  bandwidth.  We  point  out  that  the  output  distributions  can 
be  locally  negative,  as  a  consequence  of  having  chosen  the  Wigner  distribution.  This  fact  has  little 
interest  since  one  could  have  chosen  a  modified  Wigner  distribution  that  guarantees  positivity  and 
still  have  very  similar  results.  These  results  demonstrate  the  two  key  points  of  the  ABC  method: 


1.  Every  output  distribution  Cq.(i,w)  is  smoothed  in  time  and  frequency  with  respect  to  the  initial 
distribution  Cs(t,ui). 

2.  Every  stage  of  the  algorithm  emphasizes  a  different  region  of  the  characteristic  function  domain 
of  the  input  distribution. 


We  mention  two  subsidiary  issues.  First  is  the  issue  of  the  starting  distribution.  Noga,  in  his 
original  work  used  the  log-spectrogranr,  but  other  starting  distributions  are  possible,  and  perhaps 
may  be  more  effective.  How  the  final  distribution  depends  on  the  initial  distribution  can  now  be 


94 


C  (t,co)=A  (t,») 


Figure  44:  Input  distribution  Cs(t,  uj).  The  Dirac  function  is  represented  by  its  numerical  version. 


C, 


>) 


Figure  45:  Output  distribution  C\{t,u)  at  the  first  stage. 
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C2(t,<o) 


Figure  46:  Output  distribution  C2(t,ui)  at  the  second  stage. 


studied  theoretically  because  one  can  explicitly  express  the  steps  in  the  determination  of  the  kernels. 
This  should  be  explored.  Secondly,  in  the  original  formulation  the  log-spectrogram  was  used  and  this 
has  a  number  of  advantages.  We  have  not  done  so  in  the  current  formulation  because  we  wanted  to 
understand  the  time  and  frequency  averaging  operation  done  by  the  algorithm  and  keep  that  issue 
separate  from  the  benefits  of  using  the  log-spectrunr. 
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